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Index 

Upper case English characters except “I” and digit numbers indicate surfaces. 
Lower case English characters indicate coordinate systems. 


c 

w 

r, Q 

G 

P 

1 

2 

I 

n 


Matrix 


[5] 

[Lab] 


[M ab ] 




tool surface (c = G, p) 
work surface (w = 1 , 2) 
tool surfaces or work surfaces 
gear tool surface 
pinion tool surface 
pinion surface 
gear surface 
first principal 
second principal 


3 by 4 symmetric augmented matrix which relates principal curvatures and 
directions for mating surfaces 

4 by 1 matrix representing homogenous coordinates of point B 

3 by 3 matrix describing the transformation of vector from the S& coordinate 
system to S a coordinate system 

4 by 4 matrix describing the transformation of coordinates from the S& coordinate 
system to S a coordinate system 

3 by 1 matrix representing components of normal vector N 



M 


Vector 


B 

B u 

B v 

*,>'n 

h J k 


N 


n 

u 

Tl v 

v tcw) 


trV 


t? (1) v 

r* j r * 


( 2 ) 




3 by 1 matrix representing components of unit normal vector n 
3 by 1 matrix representing components of angular velocity vector u> 


position vector of point Bona surface 

dB/du 

dB/dv 

unit vectors along the principal directions of the surface at the contact point 

base vectors along axes X y 7, and Z y respectively 

normal vector of point B on a surface 

unit normal vector of point B on a surface 

dn/du 

dn/dv 

slide velocity of surfaces and 
transfer velocity 

velocity vectors of contact point in its motion over the pinion and gear surfaces, 

respectively 

angular velocity 

relative angular velocity of surface T with respect to surface Q 
tangent vector 


English Upper Case 


A 


mean pitch cone distance 


V! 


^0 j ^1 5 ^2 

.A, S 

B 

C n 

E, F, G 

Em 

E 3 

T 

I 


II 

I 


L 

L, Af, N 


M 

N 


P 

R 

Re*! E Ct 


coefficient of a quadratic equation 
auxiliary parameters 
point on a surface 
class of a function 

auxiliary parameters for first fundamental form 

machining offset 

three-dimensional space 

zero function 

first fundamental form 

second fundamental form 

Interval 

generating planar curve for a sphere 

auxiliary parameters for second fundamental form 

vector sum of machine center to back and sliding base 

middle point on the gear surface 

number of teeth 

plane 

radius of a circle 

x and z coordinates, respectively, of the center of a circle in the S c coordinate 


system 

5 coordinate system 

T the smaller absolute value of A and B 

C) the projection of V ( } on the e Cj 

Vc™ C) the projection of V iW 1 on the e Cji 

W point width 

r ^ 2 (1) the projections of vector r C (1) on vectors e 2j and e 2n , respectively 



Xmcb 

X SB 


machine center to back 


sliding base 


English Lower Case 


a 

a 


b 

b 

6 1> K 

c 


C l] > C 12 5 C 13 

d G 

c/j } 




m 


m 


FQ 

t 


S 

t 

^1 ) > t 4 

u 


u. 


constant (Chapter 1) 

semimajor axis of the contact ellipse (Chapter 3 and Appendix A) 
element of matrix [A] (i = 1,2,3 j - 1,2,3) 
constant (Chapter 1) 

semimmor axis of the contact ellipse (Chapter 3 and Appendix A) 

auxiliary variables 

clearance 

auxiliary variables 

average diameter of gear cutter 

auxiliary variables 

auxiliary variables 

gear ratio 

derivative of gear ratio with respective to <j) Q 

cradle angle 

tip radius of the cutter 

radial setting 

semimajor axis of the contact ellipse 
auxiliary variables 

surface coordinates of a cone surface 
auxiliary variables 


viii 


auxiliary variables 


U 22 > > U 32 

Non-English Upper Case 

S surface 

F shaft angle 

T angle measured counterclockwise from the root to the tangent of the path on the 

gear surface 

K ratio constant 

Sff open rectangle 

A discriminant of an equation 


Non-English Lower Case 


0 

6 


€ 




K r 


orientation angle of ellipse 
mean spiral angle 
dedendum angle 
specified tolerance value 
root angle 
principal curvature 

K + K 
/ IJ 

K — K 
I IJ 

normal curvature 

relative normal curvature of the mating surface 


ix 



A 

v i 

UJ 


<t> 

4 > 


C 

w 


$ 

■ U/ 


m\) 

*'M) 
&%(*[) 
(A^) (1) 
(A^) (2) 
V<, A^ 

T 

s 

zu 


surface coordinate of a surface of revolution 
pitch angle 

angles formed between vectors r V {1) and e 2j , and r K <2) and e 2/) respectively 
angular velocity 

turn angle of the cradle when the work is being cut 
rotation angle of the work while it is being cut 

rotation angle of one member while it is being in meshing with another member of 
a pair of gears 

transmission function, the rotation angle of the gear in terms of that of the pinion 

in a pair of meshing gears 

transmission function of a pair conjugate gear 

transmission error function 

predesigned parabolic function of transmission errors 
linear function of transmission errors induced by misalignment 
expressions of <j >\ J and A<f> f 2 in a new coordinate system 
blade angle 

surface coordinate of a cone surface and a surface of revolution 
angle measured counterclockwise from el to e 
auxiliary variable, 9 q ± <j> 
elastic approach 

angle formed by the tangent to the curvature and first principal curvature 


x 



SUMMARY 


A new approach for determination of machine-tool settings for spiral bevel gears is proposed. 
The proposed settings provide a predesigned parabolic function of transmission errors and the 
desired location and orientation of the bearing contact. The predesigned parabolic function of 
transmission errors is able to absorb piece-wise linear functions of transmission errors that are 
caused by the gear misalignment and reduce the gear noise. The gears are face-milled by head 
cutters with conical surfaces or surfaces of revolution. 

A computer program for simulation of meshing, bearing contact and determination of trans- 
mission errors for misaligned gear has been developed. 




CHAPTER 1 


INTRODUCTION 


1.1 Introduction 

The most important criteria of quality of meshing and contact of gears are the low level of noise 
and the sufficient dimensions and location of the bearing contact. Sometimes these requirements 
are contradictory and can be achieved by a compromise in the process of gear synthesis. Such a 
method of synthesis for spiral bevel gears has been developed in this report. 

Traditionally, Gleason’s spiral bevel gears are designed and manufactured with non-conjugate 
tooth surfaces. By varying machine-tool settings the transmission errors can be of different forms, 
which included a piece-wise linear function, an “S” curve, and a parabolic function, symmetrical 
or otherwise. Only a parabolic function with gear lagging is prefered. The problem encountered is 
that it is very difficult to reduce the level of a parabolic function of transmission errors with gear 
lagging. 

Litvin et al. [1] proposed a method for generation of spiral bevel gears with conjugate tooth 
surfaces. Ideally such conjugate pair provides zero transmission errors. In practice, spiral bevel 
gears are frequently required to operate under misalignment caused by mounting tolerances and 
deflections. Using the Tooth Contact Analysis (TCA) programs we have found that the conjugate 
spiral bevel gears cause lead functions of transmission errors — strong monotonous increasing or 
decreasing functions for a cycle of meshing. These functions may be considered as linear functions 
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or almost linear functions (Figure 1). Due to gear misalignment the bearing contact can be shifted 
from the desired location even to the tooth edge. For this reason it is necessary to control also the 
location and dimensions of the bearing contact. 

There is an opportunity to reach these goals if the gears will be designed as non-conjugate pairs 
that transform rotation with a predesigned parabolic function of transmission errors. Then, as it 
will be proven in the next section, a linear function of transmission errors will be absorbed and the 
sensitivity of the gears to misalignment will be reduced. 

The determination of pinion machine-tool settings is based on the local synthesis of the gears 
proposed by Litvin [2, 3, 4]. The local synthesis must satisfy the following requirements: 

1. The gear surfaces are in tangency at the chosen mean contact point. 

2. The tangent to the path of contact has the prescribed direction at the mean contact point. 

3. The contact ellipse for the tooth surfaces has the desired dimensions at the mean contact 
point. 

4. The transmission function <fr 2 {<fi 1 ) has the prescribed value at the mean contact point and its 
second derivative is negative on gear convex side and positive on gear concave side. Here, <f> 1 
and <f) 2 are the rotation angles of the pinion and gear while they are being cut, respectively. 

Requirement 4 means that the function of transmission errors is a parabolic one with gear lagging 
within the neighborhood of the mean contact point. 

Traditionally, a pair of Gleason’s gears is generated by two cones. In some cases the pinion is 
generated by a surface of revolution instead of a cone surface to obtain better bearing contact and 
to avoid an edge contact. Both cases are investigated and the machine-tool settings are determined 
according to the local synthesis and predesigned function of transmission errors. 
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1.2 Transmission Errors And Its Compensation 


In theory a pair of mating gears transforms rotation with a constant gear ratio 


m 


21 


*2 = W L 

N 2 


( 1 . 1 ) 


where w, and u> 2 are the angular velocities of the gears 

N 1 and N 2 are the numbers of teeth of pinion and gear, respectively 

Therefore, the transmission function is expected to be linear for ideal gears, i.e., 


( 1 . 2 ) 


However, the actual function <t>' 2 {<$>[) is always different from 4>' 2 {4>\) except at the mean contact 
point. The transmission errors are defined as the difference of theoretical and actual functions of 
transmission functions, i.e., 


N, 


A 4>'M) = 4>'M ) - 4>'M) = 


(1.3) 


In general the transmission errors of gears may occur due to the following four reasons [5]: 

1. The gears cannot exactly transform rotation described by equation (1.2) because of the 
method of their generation. Spiral bevel gears and hypoid gears that are generated by Gleason 
methods are good examples for this case. 
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2. The gear axes are misaligned or the gear shafts are deflected. Zhang, in his dissertation [5], 
has proved that the deflected gear shafts can be modeled as misaligned gear axes. Spur gears, 
helical gears, and conjugate spiral bevel gears are very sensitive to misalignment. 

3. Heat treatment deviation of the real gear surface is one of the most important factors in 
surface distortion. 

4. The elastic deformation of gear tooth surfaces under applied load. 

Cases 1 and 2 among the above-mentioned are the main sources of transmission errors. They 
will be discussed later. The topics of 3 and 4 are beyond the scope of this report and will not be 
discussed. 

For a pair of conjugate gears under misalignment, the investigation results in that the trans- 
mission function 4>' 3 {<t>[) becomes a discontinuous piece-wise function that is linear or almost linear 
for each cycle of meshing as shown in Figure 2. The corresponding transmission errors determined 
by equation (1.3) are also an approximately piece- wise linear function as shown in Figure 3. Such 
functions cause a discontinuity in the regular tooth meshing and usually impact at the transfer 
point. 

There is another type of function of transmission errors that is a piece- wise parabolic function 
as shown in Figure 4. This type of transmission errors does not cause a discontinuity of regular 
tooth meshing at transfer points. Gears with this type of transmission errors are not so sensitive 
to misalignment. This statement is based on an investigation into the interaction of a parabolic 
function with a linear function. 

Consider that a pair of gears is predesigned with a parabolic function of transmission errors. 
This function may be represented by 



( 14 ) 
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2 

The level of transmission errors is a(2ir/N 1 ) 

Misalignment induces a linear function of transmission errors. It may be represented by 


{^T = K 


(1.5) 


Since (A0') (l) and (A <f>' 2 ) W are very small, the principle of superposition can be applied for the 
interaction of functions (A ^) (1) and (A^'/’. Therefore, the resulting function is 

A 0; = (Ad/)"’ + (M'f' = a(<P[f + ( 16 ) 


Equation (1.6) can be rewritten in a new coordinate system by (Figure 5) 



(1.7) 


where 


At// = A 4 >' + 


4 a 



(1.8) 


From equation (1.7) we know that although the misalignment occurs, the resulting function of 
transmission errors represents the same parabolic function that has been translated with respect 
to the given parabolic function. This means that the predesigned parabolic function (A will 
absorb the linear function (A 4>[f ] induced by misalignment. The level of transmission errors 
remains the same since the parabolic function of each tooth translates the same amount. 
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Misalignment changes the path of contact. The locations of transfer points are shifted to an 
edge. The amount of the shift is determined by b/2a. In general, the absolute value of b increases 
if the amount of misalignment increases. It is possible that an unfavorable ratio b/2a may cause 
one of the transfer points to be off the tooth surface and that the function of transmission errors, 
AV’j, will become a discontinuous function for every cycle of meshing (Figure 6). To avoid this, 
the level of predesigned function of transmission error, or the absolute value of a, should be chosen 
with the expected level of transmission errors caused by misalignment. 
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CHAPTER 2 

GLEASON’S SPIRAL BEVEL GEARS 


2.1 Gleason System 

The Gleason Works, Rochester, New York, is one of the leading companies that produces 
equipment for manufacture of bevel and hypoid gears. William Gleason built the first machine in 
1874 to cut bevel gears with straight teeth [6]. During the following years, the Gleason Works has 
developed a set of machines to generate spiral bevel gears. The basic construction (Figure 7) of a 
cutting machine consists of three major parts: the frame, the cradle, and the sliding base [7, 8]. 

When cutting starts, the work is plunged into the cutter. As the cutter rotates through the 
blank, a relative rolling motion is produced between the cradle and the work spindle to generate the 
tooth surface. While the cutter rolls out of engagement with the work, the cradle reverses rapidly, 
the sliding base on which the work is mounted is translated with respect to the cutter, and the 
work is indexed ahead for cutting the next tooth. This sequence of operations is repeated until the 
last tooth is cut. 

In the process of cutting, the head-cutter rotates about its axis, and the axis generates in the 
cradle coordinate system a cylindrical surface. We may imagine that the cutter generates a tooth 
of crown gear as shown in Figure 8. Therefore, the cutting process corresponds to the motion of 
the gear rolling on a crown rack. The angular velocity of the head-cutter about its axis is not 
related with the generating motions and depends only on the desired velocity of cutting. This is 
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Top View 




Figure 8. Cutting spiral gear teeth on the basic crown rack. 
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an important advantage of the Gleason methods of manufacture. Another advantage is that the 
same method for generation can be used as well for grinding. Grinding is essential for producing 
gears with hardened tooth surfaces and of high quality. 

2.2 Head Cutters 

Traditionally straight-sided blades have been applied in practice. The blades of the cutter 
generate cone surfaces while the cutter rotates its axis. Figure 9 shows these two cones. A current 
point B on the cone surface is represented in the coordinate system S c as follows: 


' B CI ' 


r cot xp — u cos xp 

B Cy 


u sin ip sin 6 

B Cc 


u sin xp cos 6 

1 


1 


where u = ~AB and 6 are the surface coordinates, r is the tip radius of the cutter, and xp is the 
blade angle. For the inside blade of the cutter, V> is an acute angle. For the outside blade of the 
cutter, xp is an obtuse angle. 

Using equations (A.5) and (2.1) (provided usin xP / 0), we obtain the equations of the unit 


normal to the cone surface. 



n c* ’ 


sin xp 

n c = 

n Cy 

= ± 

cos xp sin 0 


_ n c z _ 


cos xp cos 9 


( 2 . 2 ) 


The total differential of vector B c is 
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- cos xp du 


[dB c ] 


sin xp [ sin 6 du + u cos d dd) 


sin xp(cos 6 du - u sin 0 dd) 


(2.3) 


The total differential of vector n c 


[dn c ] — ± 


0 

cos xp cos d dd 
— cos xp sin d dd 


(2.4) 


Equations (A. 26), (2.3), and (2.4) yield 


± cos xp cos 6 < 


— cos xp du sin xp(sin 6 du + u cos 6 dd) 


T cos xp sin 6 dd 

sin xp(cos d du — u sin 6 dd) Kim ^ ^ 


Equation (2.5) is satisfied if 


dudd — 0 


(2.6) 


One of the principal directions corresponds to du = 0; the other one to dO = 0. They can be 
represented by equations 


dB c 

de 


8B r 


d6 


e„ - 


dBc 

_ du 


8B r 


du 


(2.7) 


( 2 . 8 ) 
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Equations (2,3) and (2.7) yield 


e\ = ± 


0 

cos 6 
— sin 6 


(2.9) 


Plugging du — 0 into equations (2.5), we have 


K 

i 


1 


= T 


u tan xp 


( 2 . 10 ) 


The sense of the principal curvature relies on the chosen direction of the normal. 
Similarly, the unit vector of the second principal direction is 


'17c 


± 


— cos xp 
sin xp sin 9 


sin x p cos 6 


( 2 . 11 ) 


The principal curvature is 


= 0 


( 2 . 12 ) 


In addition to the cone surface, a tool provided by a surface of revolution is considered here. 
This surface of revolution is generated by an circular arc that rotates about the cutter axis. Such 
a surface can be applied as a grinding wheel or as a surface of a tool with curved blades. 

Suppose the generating planar curve L (Figure 10) is an arc of a circle of radius R centered 
at point (R Ct ,0,R Cx )- The spherical surface is generated by the circle in the rotational motion 
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Figure 10: Generating arc circle for the curved edge of the head cutter. 
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about the Z c -axis. Consider an auxiliary coordinate system S c > which is rigidly connected to the 
generating circle. Initially S c > and S c coincide. The generating curve may be represented in the 
coordinate system S c » with the matrix equation 


' B c . m ' 


Rc x + R cos A 



0 

B c < 

C 2 


R Cz + R sin A 

1 


1 


(2.13) 


where A is the varied parameter for planar curve L. The parameter A lies within the following 
intervals: 


Inside blade 


Outside blade 


0 < A < 7r/2 t if L is concave down; 

7r < A < 3 tt / 2, if L is concave up; 

' 3tt/ 2 < A < 27 r, if L is concave down; 
, tt/ 2 < A < 7r T if L is concave up. 


The auxiliary coordinate system S c > rotates about the Z c axis and the coordinate transformation 
from S c t to S c is (Figure 11) 


B r = 


' B Ct ' 




Be. 


— - 

H-* 

1 



10 0 0 

0 cos 6 sin# 0 

0 - sin 6 cos 6 0 

0 0 0 1 



' B c . m ' 


Be* 


Be. 


1 


(2.14) 


Equations (2.13) and (2.14) yield 
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' B Ct ‘ 


R Ct + R cos A 

Bc y 


(R c> + Esin A) sin# 

B Cl 


(E Cj + R sin A) cos 0 

1 


1 


(2.15) 


Using equations (A. 5) and (2.15), the unit normal to this spherical surface may be represented by 



' n Cx 


cos A 


n c - 

n Cy 

= ± 

sin A sin 

e 


m n Cz 


sin A cos 

e 


According to Rodrigues’ formula, the principal directions on the generating surface correspond to 
d\ = 0 and d6 = 0, respectively. The unit vector of the principal direction corresponding to dX = 0 


is 


0 


e “/ c = ± 


cos 8 
— sin 6 


(2.17) 


The principal curvature is 


sin A 

1 R Cz “I - R sin A 


(2.18) 


The unit vector of the principal direction corresponding to d8 = 0 is 


± 


— sin A 
cos A sin 6 
cos A cos 9 


(2.19) 
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The principal curvature is 


K 


a 



( 2 . 20 ) 


2.3 Coordinate Systems and Sign Convention s 

Left-hand gear-members are usually cut by the counterclockwise motion of the cradle that 
carries the head-cutter. This motion is viewed from the front of the cradle and from the back of 
the work spindle. Cutting is performed from the toe to the heel. Figure 12 shows the top and front 
views of the machine when a left-hand gear-member is cut. 

Right-hand gear-members are usually cut by motions that are opposite to the motions of the 
left-hand members being cut. Cutting is performed from the heel to the toe. Figure 13 shows the 
top and front views of the machine for this case. 

We set up five coordinate systems in either case. Coordinate system S c is rigidly connected to 
the head cutter, coordinate system S w is rigidly connected to the work, and coordinate systems 
Sm, S p and S a are rigidly connected to the frame. Axes Z m and Z p coincide with the root line and 
pitch line, respectively. Axis X m is perpendicular to the generatrix of the root cone of the work. 
Axis X p is perpendicular to the generatrix of the pitch cone of the work. Axes Z a and Z w coincide. 
Origin Om is located at the machine center, and origins O a and O v are located at the apex of the 
pitch cone of the work. 

Three special machine-tool settings, which are the machining offset, machine center to back, 
and the sliding base, are used only for the generation of pinions. The machining offset, denoted 
by E m , is the shortest distance between the cradle axis and pinion axis. In figures 12 and 13, T m 
represents a vector sum of machine center to back, X M cb, and the sliding base, X SB . The change 
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TABLE 1: SIGN CONVENTIONS OF MACHINE-TOOL SETTINGS. 



Right-Hand Member 

Left-Hand Member 

Cradle Angle 
Q 

+ 

counterclockwise (CCW) 

clockwise (CW) 

- 

clockwise (CW) 

counterclockwise (CCW) 

Machining Offset 

E rn 

+ 

above machine center 

below r machine center 

- 

below machine center 

above machine center 

Machine Center to Back 

Xmcb 

+ 

work withdrawal 

work withdrawal 

- 

work advance 

work advance 

Sliding Base 

X SB 

+ 1 

work withdrawal 

work withdrawal 

- 

work advance 

w r ork advance 

Lrn 

+ 

A 'sb - + and Xmcb '■ ~ 

Xsb '■ + and Xmcb - - 

- 

X sb • ~ and Xmcb • + 

Xsb : “ an d Xmcb : + 


of machine center to back is directed parallel to the pinion axis and the direction of the sliding base 
is pointed parallel to the cradle axis. 

The sign conventions for machine-tool settings are given in Table 1. 

2.4 Generated Tooth Surfaces 

The generated surface Ew is an envelope of the family of the tool surface Ec- Surfaces Ew 
and Ec contact each other at every instant along a line which is a spatial curve. Surface E^ is 
conjugate with E c- In mathematical sense the determination of a conjugate surface is based on 
the theory of an envelope of a family of given surfaces. In differential geometry, to determine E w 
we must find: 

(a) the family of surfaces E$ generated by the given surface E c in the S w coordinate system 
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and 


(b) the envelope Ew of the family of surfaces E$. 

The matrix representation of the family of surfaces £$ may be represented by the matrix equation 


\B W ) = [M wc ] [ B c ] 


(2.21) 


where [M wc ] is a matrix which describes the transformation of coordinates from the “old” coordinate 
system S c to the “new” coordinate system S w . From Figures 12 and 13, we obtain 


[M wc ] = [M wa ] [M ap ] [Mpm] [Mmc] 


We can obtain [M wa ] from Figure 14 as 


[M wa ] = 


cos <f> w ± sin (f> w 0 0 

T sin <f> w cos <f> w 0 0 

0 0 10 

0 0 0 1 


(2.22) 


(2.23) 


where (f> w is the rotation angle of the work while it is being cut. Here the upper sign corresponds 
to the generation of a left-hand spiral bevel gear that is shown in Figure 12, and the lower sign 
corresponds to the generation of a right-hand spiral bevel gear shown in Figure 13. Henceforth we 
will obey this notation. 

The transformation matrices [M ap ] and [Mpm] can be obtained from Figures 12 and 13 as 
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cos fi 0 sin fi 0 


[Map] = 


0 10 0 
— sin // 0 cos fi 0 


[M pm ] = 


[ 0 0 0 1 J 

cos 8 0 — sin 8 — L m sin 8 

0 1 0 ±E m 

sin 6 0 cos 8 L m cos 8 

0 0 0 1 


(2.24) 


(2.25) 


where ji and 8 are the pitch angle and dedendum angle of the work, respectively. To derive the 
transformation matrix [M mc ], let us apply an auxiliary coordinate system S 8 rigidly connected to 
the tool (Figure 15). Thus 


[Mmc] = [M ma ][M 3C ] 


' 1 

0 

0 

0 ' 


' 1 

0 

0 

0 

0 

cos <f> c 

± sin <f> c 

0 


0 

cos q 

^ sin q 

s sin q 

0 

T sin 4> c 

cos 4* c 

0 


0 

±sin <7 

cos q 

s cos q 

0 

0 

0 

1 


0 

0 

0 

1 


where <f) c is the turn angle of the cradle while the work is cut, and s is the radial setting. The 
determination of the envelope Tiw of the locus of surfaces £$ is based on necessary and sufficient 
conditions of envelope existence that have been developed in the classical Differential Geometry. A 
simpler method representation for determination of necessary condition of Yjw existence is based 
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on the geometric property of conjugate surfaces: at points of tangency of the generating surface 
E c and the generated surface the common unit normal n to the surfaces is perpendicular to 
the slide velocity G f these two surfaces [9, 10]. This is given by the scalar product 

n-F <CW) = 0 (2.27) 

In the modern theory of gearing, equation (2.27) is called the equation of meshing. This equation 
is of fundamental importance in the kinematics of gearing. Since equation (2.27) is valid in any 
reference system, we will derive the equation of meshing in the S m coordinate system. Let us 
designate t r V^ } and tr V , ^ the transfer velocities of common contact points B m on the cutter and 
the work, respectively. Thus 


- t7 y 

Y m tr r 


(C) 

m 


trV, 


(W) 


(2.28) 


The cradle rotates about the Z m axis with the angular velocity (Figures 12 and 13); therefore, 
the transfer velocity tr V m is represented by the equation 

trV^ = c£ x B m (2.29) 


The work rotates about the Z a axis with the angular velocity u5 m (Figures 12 and 13) which does 
not pass through the origin O m of the £ m coordinate system. It is known from the theoretical 
mechanics that the angular velocity w rn may be substituted by an equal vector u; m which passes 
through Om and a vector-moment 

(2.30) 


OmOn X UJ. 
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Note that the moment has the same unit and physical meaning as linear velocity. Thus 


tr 


v iW) 


-(WO 


x B rn 4“ O m O a X w 


(W) 


(2.31) 


It is evident from Figures 12 and 13 that 


1 



(C) 




0 


0 


(2.32) 



— sin 7 

0 


cos 7 


O m O a = 


0 

TEt, 


(2.33) 


(2.34) 


In equation (2.33) 7 is the root angle of the work. Substituting equation (2.29)-(2.34) into equation 
(2.28), we obtain 


V lcw) 

r m 


u> [W) (E ' m ±S m ,)cos 7 

±u> {C> B mz +u’ ,W) [(5 mi T L m ) sin 7 ^ B mt cos 7] 
±w (W) (5 my -£ m )sm 7 


(2.35) 
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The coordinates of the common contact points i? m may be obtained from equations of the generating 
surface Sc- Then we get 


[B m ] = [Mmc\ [B e ] 


(2.36) 


The common unit normals n m may be represented by the unit normals to Sc- Therefore 


[tt-m] — [fmc] [n c ] 


(2.37) 


where [L mc ] is the rotation matrix obtained by eliminating of the last row and last column of the 
corresponding matrix [M mc ]. 

Hence, if Sc is a cone surface, substituting equations (2.1) and (2.26) into equation (2.36) we 
obtain 


' B mt ' 


r cot xp — u cos ip 

B my 


u sin ip sin r =f= s sin (q - <p c ) 

B mt 


u sin xp cos r -f s cos($ — <p c ) 

1 


1 


where r = 6^q±<p c . Substituting equations (2.2) and (2.26) into equation (2.37) the unit normals 
may be represented as 


^m z 


sin xp 

^m y 

= ± 

cos xp sinr 

_ 


cos xp cos r 


(2.39) 


Similarly if Sc is a spherical surface, substituting equations (2.15) and (2.26) into equation (2.36), 
we obtain 
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r 

w 

e 

cq 


R Cz + R cos A 

Bm y 


( R Cl — R sin A) sin t ^ s sin(g - <f> c ) 

B mz 


{R Ct + R sin A) cos r + s cos (q - (p c ) 

1 


1 


Substituting equations (2.16) and (2.26) into equation (2.37) the unit normals may be represented 
as 


n m x 


cos A 

n TTLy 

= ± 

sin A sin r 

z \ 


sin A cos r 


(2.41) 


Designate 


(C) 




U ) 


(W) 


(2.42) 


Using equations (2.27), (2.35), (2.38), and (2.39), we may obtain the equation of meshing for the 
case that £c is a cone surface by 


[u — r cot xp cos xp) cos 7 sin r + s[(m cvr — sin 7 ) cos xp sin 6 =p cos 7 sin xp sin(g — <p c )] 
±i? m (cos 7 sin xp + sin 7 cos xp cos r) — Z m sin 7 cos xp sin r - 0 


(2.43) 


For Yjq being a spherical surface, using equations (2.27), (2.35), (2.40), and (2.41), the equation of 
meshing is represented by 


( R Cz cos A - R Cx sin A) cos 7 sinr + s[(m cvr - sin 7 ) sin A sin 0 4 = cos 7 cos A sin(</ — <p c )] 
±£ m (cos 7 cos A + sin 7 sin A cos r) - L m sin 7 sin A sin r = 0 


(2.44) 
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Equations (2.43) and (2.44) relate the generating surface coordinates (u and 6 for a cone surface 
or A and 9 for a surface of revolution) with the turn angle 4 > c . 
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CHAPTER 3 


SYNTHESIS OF SPIRAL BEVEL GEARS 


3.1 Gear Machine-Tool Settings 

We designate for the following discussions the gear-generating tool surface by £ C , the generated 
gear surface by £„ the pinion-generating tool surface by Ep, and the generated pinion surface Si- 
A parameter with the subscript i indicates that it is related to surface E„ To set up the gear 
machine-tool settings, the following data are considered as given: 


T: shaft angle 

jV 2 : gear tooth number 

JV ; pinion tooth number 

^ : gear root angle 

A: mean pitch cone distance 

f3\ mean spiral angle 

xp G : blade angle for gear cutter 

d • average diameter of gear cutter 

W G : point width 
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3-1.1 Preliminary Considerations 


We prefer to calculate the valuer of pitch angles and dedendum angles rather than obtain them 

from the blank design summary because the data in the summary are not accurate enough for 
computer calculations. 

The gear pitch angle is represented by 


The pinion pitch angle is 


/i 2 = arc tan 


sinT 


N, 

N. ' 


~ + cosT 


mi = r - m 2 


(3.1) 


(3.2) 


The dedendum angles are 


- Mi - ?j , S 3 = fj , 2 - 


(3.3) 


3-1.2 Gear Cutting Ratio 

The process of gear generation is based on the imaginery meshing of a crown gear with the 
member-gear. The instantaneous axis of rotation by such meshing coincides with the pitch line, 
axis that is shown in Figures 12 and 13. The generating surface E c , which may be imagined as 
the surface of the crown gear, and the to be generated gear surface St contact each other at a line 
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at every instant. The ratio of angular velocities of the crown gear and the being generated gear (the 
cutting ratio) remains constant while the spatial line of contact moves over surfaces S G and S 2 . 
The determination of cutting ratio is based on following consideration. The angular velocity in 

relative motion is 


G2) ^(G) ^(2) __ 7 

(jj — ijJ — U/ — drop 


(3.4) 


This means that vectors J l ° 2) and k p are colhnear. Since equation (3.4) is valid in any reference 
frame, let us derive it in the S m coordinate system. Prom Figures 12 and 13 we have 


*p = 


sin S 2 

0 

cos 6, 


(3.5) 


By replacing the superscript ‘c< b, 'o’ in equation (2.32) and •»' by '2' in equation (2.33), we may 
represent in matrix from angular velocities and Consequently, we obtain the following 

equation 

{ n\ 1 21 ( 2 ) 

(3.6) 


q:ui (C| f/ sin 7 Z _ W* cos 7 2 


sin 


cos < 


Equation (3.6) results in that 


m 


G2 


(G) 

u> sm 

o; (2) COS S 2 


(3.7) 
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3 ' 1 - 3 Cutter Tip Radius, Radial Setting, and Cradle Angle 


Figure 16 shows that the inside and outside tip radii of the head- 


cutter are represented by 


r o = ~{d G T w 0 ) 


(3.8) 


Figure 16 show, the front view of the installation of the head cutter. From the relation, between 
the lengths and angles of the triangle 0 m 0 c M o , we may express the radial setting s c and cradle 


angle q G as follows: 




G ~\ 4 + ^ cos $2 - d G A cos S 2 sin f3 


(3.9) 


and 




arccos 


A 2 2 r 2 

A cos 6. + S i 

2 G 4 


2 As g cos S 2 


(3.10) 


3.2 D etermination of the Mean Contact Pgint o n the Gear Tboth Surfaces 

The gear and pinion surfaces of spiral bevel gears are in point contact at ever, instant. The 
mean contact point is the center of the bearing contact and its location is selected generally at the 
middle of the working depth on the gear tooth. Figure 17 shows a gear tooth surface. Section AD 
■s the gear tip and it i, parallel to the generatrix of the root cone of the pinion. Section W' is the 
pinion tip and it is parallel to the root line of the gear. The working area is within D ABCD. In 
the S p coordinate system, line ~AD may be represented by 
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Figure 16: 



Figure 17: Gear tooth surface. 
42 


B Vt = B Pz tan 8 1 - c 


(3.11) 


where c is the clearance and is the pinion dedendum angle. Line BC is represented by 

B Pt = -B Vl tan + c (312) 


The mean contact point is located on a line which passes through the middle point of the two points 
at which the normal section of the gear surface intersects line AD and line BC, respectively. In 
addition, the mean contact point must be on the gear surface. This means that it must satisfy the 
equation of meshing for the gear being generated by the tool. We will use these two requirements 
to determine the location of the mean contact point and represent the procedure of derivations as 
follows 

Step 1: The initial guess for 6 C is 


6 g = ±{q G - 0 + * 72 ) 

Step 2: Determination of u G based on the given 0 G 

Equation (2.43) determines parameter u G . The turn angle 4> a is set to zero when equa- 
tion (2.43) is applied. 

Step 3: Representation of gear tooth surface in coordinate systems S c and S p 

Equation (2.1) determines the gear tooth surface in coordinate system S c . The gear tooth 
surface may be represented in the S p coordinate system as follows: 


[B p ] = [Mpm] [ Mmc } [Bo] 


(3.13) 
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Transformation matrix [Mpm] may be obtained from equation (2.25) by setting E rn and L m 
to zero. Equation (2.26) determines matrix [M mc \. 

Step 4: Determination of middle point 

The x coordinate of the middle point M of the two points, which are the intersections of the 
normal section of the gear tooth surface and the gear tooth tips, may be obtained by 

Ur . - -«■“«.) (3 . 14) 

The above equation is derived by dividing the sum of equations (3.11) and (3.12) by 2. 

Step 5: Judgement of u G 

The acceptable value of u G is determined by the following criterion: 

\ B Pr - I < f 


where e is a specified tolerance value. If the above criterion is satisfied, parameters u G and 0 G 
of the mean contact point are determined. Otherwise, repeat Step 2 to Step 5 by a new 
value of 0 G until the criterion is satisfied. 

As a matter of fact, the determination of the location of the mean contact point is the same as 
that of a root of equation 


B Vv — M Vz — 0 


(3.15) 
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The new value of the 6 G in Step 5 depends on which method is used to solve this equation. In 
this study Newton's method was used. 

So far we have already determined parameters u G and 6 G of the mean contact point. Repeating 
the task done in Step 3, we have the coordinates of the mean contact point B. The common unit 
normal n to surfaces E G and E 2 at the mean contact point B is 


[ n p] — [^pm] [^mc] [«c] 


(3.16) 


where matrix [I pm ] is obtained by deleting the fourth row and column from matrix [MpJ given 
by equation (2.25). Similarly, we may obtain rotation matrix [L mc \ from matrix [M mc \ by equa- 
tion (2.26). Although the unit normal has two directions, we choose the direction corresponding to 
the positive sign in equation (2.2) regardless of the hand of the gear. The principal directions at 
the mean contact point B on the gear tool surface E G in the S p coordinate system may be obtained 
by the following coordinate transformation: 


'ui v 


[Lprn] [fmc] e 


(3-17) 


Here we choose positive sign in equation (2.9) as the direction of the first principal direction. 
The second principal direction is determined by rotating of the first principal direction about unit 

normal by 90°. 

The principal curvatures and directions at the mean contact point B on the gear surface E 2 may 
be derived according to the formula expressed in Section A.2. Note that surfaces E G and E 2 are 
in line contact. To apply these formula, we may consider that surfaces E 2 and E G are equivalent 
to surfaces E^ and Eg, respectively, in Section A. 2. 
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The derivation of the principal curvatures and directions at the 
formed as follows: 


mean contact point 


Step 1: We represent the angular velocity u> W in the S p coordinate system as follows: 


< 2)1 , ( 2 ) 


- sm fi 2 

0 

[_ cos fl 2 


This is a direct result from drawings of Figures 12 and 13. 

Step 2: We represent the angular velocity Cj G) in the S p coordinate system as follows 


(0)1 __ <G) 


= 


cos 6 , 


sin£„ 


This is also a direct result from Figures 12 and 13. 
Step 3: The relative angular velocity w* a) is represented by 


,-( J 0) -#(2) _(G) (2) 

U> p = U) p - Ui p = ±U> 


— sin fi 2 + m G2 cos 8 2 


cos fi 2 + m G2 sin S 2 


Step 4: The transfer velocity of the mean point B on surface is 


T?(0) -(O) S 

tr^ p = W p X B p 


B is per- 


(3.18) 


(3.19) 


(3.20) 


(3.21) 
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Step 5: The transfer velocity of the mean point B on surface E 2 is 


Step 6 


Step 7 


Step 


step 9 


t 7< 2 » -< 2 > v D 

tA' = 


(3.22) 


: The relative velocity of the mean point B is 


t --( 2G > 




(3.23) 


the projection of V^ G) on the e 


is 




(2G) t 7(2G) - 

/ =V P 


(3.24) 


The projection of V p on the 


is 


T/ (>o) _ t 7( 2G ) - 

^ g /7 —^p ' e <3j T r 


(3.25) 


: Using equation (A. 33), we obtain 


\A 2a) 

’ K Gj V °1 


U/' 


.(2G) . 


U P e G, 


(3.26) 
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Step 10: Using equation (A. 35), we have 


ui 2G) - 

t v G n 


-*( 2 G) _ 


(3.27) 


Step 11: Using equation (A. 36), we obtain 




- J2G) 

n pU> p 


V (2G) 

p 


Tin 




.(2) 


x trV 


(G) 


JG) 


r“r( 2) 1 


trV 


(3.28) 


Note that m f =0 

G2 

Step 12: To determine the principal directions at point B on gear surface, we first use equa- 
tion (A. 40). Thus 


tan2cr 2G = 


2 fl 13 fl 23 


a 23 - a n +(* G/ - * G/i K 


(3.29) 


Rotating unit vector e Gj about the unit normal vector n by -cr 2G , we may obtain unit 
vector € 2 j . Rotating unit vector e 2 ^ about the unit normal vector ti by tt/2, we may have 
unit vector el . 

2 II 

Step 13: Using equations (A. 41) and (A. 42), we may determine the principal curvatures on the 
gear surfaces as follows: 


= 


Q 23 ‘ fl 13 


+ {k 


G I 


,K 


a 33 COS 2£r 2G 


(3.30) 


K 2j + — ( K Gj + K a n ) 


2 2 
a + a 

1 O *» 


(3.31) 
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Step 14: Eliminating k 2 by considering the sum of equations (3.30) and (3.31) together and then 
dividing the sum by 2, we can determine Eliminating by dividing the difference of 
equations (3.31) and (3.30) by 2, we can determine k 2 . 

3.3 Local Synthesis 

The determination of pinion machine-tool settings is based on the idea of local synthesis of 
gear tooth surfaces proposed by Litvin ;2, 3, 4]. The goal of local synthesis for meshing of spiral 
bevel gears is to satisfy the following requirements: 

1. The gear tooth surfaces must contact each other at the prescribed mean contact point B. 

2. The contact ellipse for the gear tooth surface must have the desired dimensions at point B . 

3. The tangent to the contact path must have the prescribed direction at point B. 

4. The instant gear ratio m 2l (^j) and its derivative m' (<£ x ) must have the prescribed values at 
point B. 

The local synthesis for the gear tooth surfaces connects the concept of meshing and the concept 
of bearing contact. It provides the optimal conditions of meshing for the gear tooth surfaces being 
in mesh at, and within the neighborhood of, the mean contact point B . The local synthesis needs 
the information on the characteristics of the tooth surfaces of the zero, first, and second orders. 

Starting the local synthesis we already know the location of the mean contact point B on the 
gear surface, the unit normal to the gear tooth surface at point #, the principal curvatures and 
directions at point B on the gear tooth surface. 

We will consider the local synthesis in a fixed coordinate system S y . Figure 18 shows the 
relations among 5/, fixed coordinate systems which are attached to the frame of the gear generator, 
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and fixed coordinate systems connected to the frame of the pinion generator. From Figure 18 we 
know that Sf and S p (G) coincide with each other. Therefore, the coordinates of the mean contact 
point B : the orientation of the surface unit normal, and the principal directions at the point B on 
the gear surface are known since they have been determined in the 5 (oj system. 


3.3.1 Preliminary Considerations 

Spiral bevel gears transform rotation motion between intersecting shafts with an instantaneous 
point contact of surfaces. It corresponds to the second case discussed in Section A. 2. Some elements 
of matrix [A] shown in equation (A. 30) are not related with the principal curvatures and directions 
of the pinion surface; therefore, they may be derived at the stage where the principal curvatures and 
directions are not known yet. We will consider that all derivations are performed in Sf coordinate 
system. Throughout the rest of the report, we will drop the subscript if it is considered in the Sf 
coordinate system. 

The following representation is the result of direct observation of drawings of Figure 18. 


mi , (i) 

u> = ±u> 


sin /ij 

0 

L cos/!* J 


(3.32) 


( 2)1 , 0 ) 
u) = ±m 01 u> 


— sin fi 2 

0 

L cos M 2 J 


(3.33) 


Recall that the upper sign in the equations corresponds to a left-hand member. As far as a pair 
of spiral bevel gears is concerned, the hands of the spiral must be opposite; a left-hand gear 
(pinion) and a right-hand pinion (gear) constitute a pair. Therefore, if we take the upper sign in 
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Left-Hand Gear 
Right-Hand Pinion 





equation (3.32), we must pick up the lower sign in equation (3.33). The relative 


-< l2 > -(1) _(2) 
LU — U) — jjJ 


The transfer velocity of the mean contact point B on surface Ej is 


trV {V> = 


w 1 ’ x B 


The transfer velocity of the mean contact point B on surface S 2 is 


T ?< 2 ) _(=) s 
t T V = u x B 


The relative velocity of the mean point B is 




r-" = tr v ll) - tr v 


;(2) 


The projection of V* * on the vector e"„ is 

2 I 


_ -( 12 ) „ 
V *I = V ■ e 2 


The projection of V ( 1 on the vector e. is 

2 n 


V (12) -( 12 ) _ 
V2 n = V • e 


2 n 


■ngular velocity is 


(3.34) 


(3.35) 


(3.36) 


(3.37) 


(3.38) 


(3.39) 
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Let surfaces Ei and Ejf, E 2 and Eg, be equivalent, respectively. Equation (A.33) yields 


(12) 




r-(l2) - 

k 


ne 




(3.40) 


Using equation (A. 35), we obtain 


( 12 ) 






-(12) — 
ne. 


(3.41) 


Equation (A. 36) yields 


/ ( 12 ) \ (, ,r( ,2 )\ r -- (12) l? (1J, l 

«2, (< ) + *2„ (V 8 ) - [«“> ^ J 


(3.42) 




J" x tr V' <2) - w*' x trF 0 ’) + (^' 1) ) m' ai (n x k 2 ) ' B 


where k 2 is the unit vector along the axis of rotation of the gear. It is represented by (Figure 18) 


— sin fji 2 


N = 


o 


COS /i 2 


(3.43) 


In general, spiral bevel gears are designed and manufactured with non-conjugate tooth surfaces. 
Varying the machine-tool settings it is possible to obtain a lead function of transmission errors, a 
parabolic function with pinion lagging, or a parabolic function with gear lagging. Only a parabolic 
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function with gear lagging is good for applications. Therefore, for the convex side of gear tooth m' 2i 
we must provide a negative value, and for the concave side of gear tooth m' 2i must be positive. The 
absolute value of m' 2i controls the level of the transmission errors. We consider m' 2i as an input. 

On the gear surface a path of contact that appears almost straight and substantially vertical 
to the root may fully satisfy the operating requirements in many cases; however, it should not be 
assumed that this is true for all cases. Sometimes a different direction or shape may be preferable 
[11]. The tendency of the direction of the contact path may be determined by the relative veloc- 
lty r V at the mean contact point on the gear surface. Let v 2 denote the angle between the unit 
vector e 2j at the mean contact point on the gear surface and the direction of tangent at the same 
point to the path of contact. The relation between the principal directions and the direction of the 
contact path may be represented as follows (Figure 19): 


v 2 = T + a 


2 G 


(3.44) 


The angle T is measured counterclockwise from the root to the tangent of the path. This angle is 
considered as an input. 

3 ' 3 - 2 Relations Between Directions of the Pat h s of the Mean Contact Point in its Mnfi on over the 
Gear and Pinion Tooth Surfaces 

Figure 20 shows the common tangent plane to the gear and pinion surfaces at the mean contact 
point B. The notations in Figure 20 are as follows: 

e 2; and e 2jj : unit vectors of the principal directions on the gear surface 

-*(12) 

V sliding velocity at point B 

■ velocity vector of contact point B in its motion over the pinion surface 

-*(2) 

rV [ velocity vector of contact point B in its motion over the gear surface 
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Figure 20: Common plane at the mean contact point. 
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r v£’ and rtf": the projections of vector r V (1) on vectors e 2; and e 2/j 

n l and u 2 \ angles formed between vectors r l /l 1 and e 2j , T V and e 2] , respectively 

The relation between angles v x and v 2 depends on parameters in motion and the principal 
curvatures of the gear tooth surface. For derivations we will use the following equations: 


r F (J> = r y (1) + y ,12) 


(3.45) 


that yields 


A 2 J — rV2 j -r V '2 J 


(3.46) 


,,(») t ,<>) , T r< 12 > 

A 2 n — A 2 n “t" V2 n 


(3.47) 


From the geometric relations shown in Figure 20 we have 


t J2) _■.(*). 

rViu = rV 2 , tan v 2 


(3.48) 


,.(» T ,(») . 
r Vi n = rVj, tan V, 


(3.49) 


Substituting equation (3.48) and (3.49) into equation (3.47), and then substituting equation (3.46) 
into (3.47), we obtain an expression for r V 2j * in terms of V 2j , V 2// , v x , and v 2 as follows 


_ • ■‘17 

»**/ - ~T. 


Vi? - V 2 \ l2) tan n 2 


tan v 2 - tan iq 


(3.50) 
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According to equation (A. 29), r V 2j and rVjjJ are related as follows: 


tr (1) , 

°31 r' 2 j + a 


,/■) _ 

32 t"2jj — fl 33 


(3.51) 


Here surface E 2 is equivalent to surface Eg; surface Ej to surface E^. Substituting equation (3.49) 
into equation (3.51), we have 

( a 3i + a 32 v l ) T V^ ] = a 33 ( 3.52) 


Finally, combining equations (3.50) and (3.52), we have the relation between angles v x and v. 


tan ^ = 


(a33 + a3 1 < 12 ’)tan^-a 31 F 2 ( i 2) 


.O' 


r(12) 
2 II 


- F'; 21 tan v 2 ) + a 3 


(3.53) 


^ ^ ^ Principal Curvatures and Directions of the Pinion Tooth Surface at the Mean Contact Point 

The derivation of principal curvatures and directions of the pinion tooth surface at the mean 
contact point is based on the following procedure. 

Step 1: Representation of A and B in terms of coefficients a u , a 12) and a 22 

We recall that the lengths of semiaxes of the contact ellipse, a and 6, are determined by 
parameters A, B, and e(see Section A.4). 

The sum of equations (A. 31) and (A. 34) yields 
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fl H + G 22 K 2S 


( 3 . 54 ) 


Substituting equation (A. 31) by equation (A. 34) we obtain 


- fC 1A cos 2 <t 12 


(3.55) 


We may represent parameter A in equation (A. 54) in terms of a n , a 12 , and a 22 as follows 


A = - 


1 

4 




- a 22 ) + 4(1 


(3.56) 


Also, the representation of parameter B in equation (A. 55) in terms of a n , a 12 , and a 22 gives 


B= - 


1 

4 




- a 22 ) + 4a 


(3.57) 


Furthermore, equations (3.56) and (3.57) yield 


[(a u + a 22 ) + 4.4] = (a n - a 22 ) + 4a 12 - [(a 


a 22 ) + 45] 


(3.58) 


Let T denote the smaller absolute value of A and B. Therefore, equation (3.57) can be written 
as 


[(a u + a 22 ) + 4T] =(a n -a 22 ) + 4a 1: 


(3.59) 
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STEP 2: Representation of coefficients a n , a l2 , and a 22 in terms of r V 2 ( ^ and r V 2 ( ^ 

Using the first two equations in (A. 29) and equation (3.54), we may derive a system of three 
linear equations in unknowns a n , a 12 , and a 22 


i tt-U) 

2 / a u + rVi r/ a, 


= a. 


T/f 1 ) I xrO) 

r*2j a i2 & 2 


= a. 


(3.60) 


a n + 


+ 


K 


A ; 


where k; a — ac 2S • Using Cramer s rule we may solve equations (3.60) as follows: 


a, i = 


- »„ K' + «. 
K’)‘ + (,<’)’ 


(3.61 ) 


T r(0 , 

a u T* 2 n + a 2 


t rO ) t r( 1 ) T /■( 1 ) 

rKj ~ rU 2/ r V 2 n 


{<) +{,<')' 


T Al) , 

-a 13 r l 2 7 + a. 


» + « A (,<•’)■ 


i-o +«’)' 


The third equation in (A. 29) is 


T/ (i) (i) 

a 31 i - * 2 / "h fl 32 — U 3 


(3.62) 


(3.63) 


(3.64) 


Substituting equation (3.49) into equation (3.64), we have 


F (1) - 

r*2j — 


a i3 + a 23 tan V 1 


(3.65) 
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Plugging equations (3.49) and (3.65) into equations (3.61)-(3.63), 
results 


Q 11 = d i k a + 


G 12 = d 2 K A + 6 2 


a !3 = d 3 K A + 


where 




2 

tan i/j 
1 + tan 2 v x 


— tan v x 
1 + tan 2 v x 


d 


3 


1 

1 + tan 2 v j 


ft, = 


<4 ~ 4 tan ' y i 

fl 33 (1 + tan ' V l) 


( a 23 + Q u tan i/, )(a 13 + a 23 tani/J 
033(1 + tan 2 iq) 


obtain the following 


(3.66) 


(3.67) 


(3.68) 


(3.69) 


(3.70) 


(3.71) 


(3.72) 


(3.73) 
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Step 3: Determination of k a 

A 

Equations (3.59) and (3.66)— (3.71) lead to 

[4T 2 -(!>;+ 6^)] (1 + tan 2 
2T(1 + tan 2 v x ) + b 1 (1 — tan 2 ) + 2 b 2 tan v x 


(3.74) 


Since k a = K 2r - k iE , equation (3.74) becomes 


4 T -(b[ + %) (1 + tan v x ) 


2T(l + tan i / l ) + <^(1 - tan i/J + 26 2 tantq 


(3.75) 


Note that 



(3.76) 


where t is the semimajor axis of the contact ellipse. This is an input datum. In general, it 
is about one sixth of the width of the gear tooth. Gleason Works suggests that the elastic 
approach e is 0.00025 inches [11]. 

Step 4: Determination of a u , a 12 , and a 22 

Substituting equation (3.74) into equations (3.66)— (3.68), we obtain a n , a 125 and <z 22 . 

Step 5: Determination of cr 12 

Using equations (A. 32) and (3.55), we obtain 


tan2cr 12 = 


^ a i2 

k 2A — a n + a 22 


(3.77) 
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It provides two solutions for a 12 , and we will choose the smaller value. Rotating unit vector e 2j 
about the unit normal vector n by — cr l2) we may obtain unit vector e, . Rotating unit 
vector e x about the unit normal vector n by tt / 2 , we may obtain unit vector e 2 ^ . 

Step 6: Determination of k ia 

Using equation (A. 32), we obtain 


sin 2 <7 12 


(3.78) 


Step 7: Determination of k, and k. 

The principal curvatures of the pinion surface at the mean contact point B are determined 

by 


= 


'j - 


K. = 
1 II 


(3.79) 


3.3.4 First Order Characteristics 

Four surfaces, the gear head-cutter surface Eg? the gear surface E 2 , the pinion head-cutter 
surface Ep, and the pinion surface Ei, are in tangency simultaneously at the mean contact point B. 
It implies that these four surfaces have a common normal at the mean contact point. We can use 
this information to determine pinion blade angle ip p and parameter r p . 

The representation of the unit normal to the pinion head-cutter surface in the 5/ coordinate 
system is 
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l n f ) - { L fp(f , )}[L p (p) rn (p))\L Tni P)^P)][n ci p ) ] 



■ -1 

0 

0 ' 


cos S l 

0 

— sin 

= 

0 

-1 

0 


0 

1 

0 


0 

0 

1 _ 


sin 8 j 

0 

cos 8 j 


r i 

0 


0 

i r 

1 0 


0 cos </> p ± sin 4> p 

0 ^fsin<^p cos <p p 


0 

0 cos q p ^sin^ F 
0 ± sin q p cos q p 



’ n c (P| " 


n JP) 

c y 

_ 

1 

' w b 

e 

- - i 


(3.80) 


Let us consider the straight-edged blade first. Equation (2.2) describes the unit normal in 
the S c coordinate system. Before plugging equation (2.2) into equation (3.80), we must investigate 
the sense of equation (2.2). From Figure 18 we know we must choose the minus sign for the unit 
normal. Therefore, equations (2.2) and (3.80) yield (subscript ‘f’ is dropped) 


" n x 


cos 8 X sin xp p — sin 8 j cos xp p cos r p 

Tly 

- 

cos rfip sin r p 

. n z . 


— sin 8 l sin xp p - cos 8 1 cos x \> p cos r p 


(3.81) 


Multiplying n x by cos£ 1? n z by — sin 8 li and then considering their sum, we obtain 


n x cos 8 Y - n 2 sin 8 1 = sin xp p 
Obviously, the pinion blade angle is 


p - < 


arcsin(n a; cos 8 X - n z sin S l ) 

(tt - V-p) 


Gear Concave Side 
Gear Convex Side 


(3.82) 


(3.83) 
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The x component in equation (3.81) may be rewritten as 


COS T p 


n x — cos 8 X sinipp 
— sin 8 X cos xp p 


The y component in equation (3.81) may be rewritten as 


sin r p = 


cos ip 


p 


The parameter r p may be obtained by 


sin t p 

r p — 2 arc tan — 

1 + COS T 


Let us now consider the curve-edged blade. Equations (2.16) and (3.80) yield 


’ n x " 


cos 6 X cos A p — sin S x sin A p cos r p 


— 

sin Ap sin r p 

- n z 


— sin 6 X cos Ap — cos 6 X sin A p cos r p 


Multiplying n x by cos^, n z by — sintfj, and then considering their sum, we obtain 


cos A p = n x cos b x — n z sin 6 X 


(3.84) 


(3.85) 


(3.86) 


(3.87) 


(3.88) 
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The quadrant in which the parameter A p locates may be determined by the discussion stated in 
Section 2.2. 

The blade angle is the angle formed by a line tangent to the blade surface at the mean contact 
point and a line perpendicular to the cutter head face. Thus we have 


f 5/2 tt - A p 


= 


3/2t r - Ap 
1/2tt - Ap 
3/2tt - Ap 


pinion concave side, blade concave down; 
pinion concave side, blade concave up; 
pinion convex side, blade concave down; 
pinion convex side, blade concave up. 


Rewriting the x component in equation (3.87), we have 


cos r p ~ 


n x — cos 6 1 cos A p 
— sin sin A p 


(3.89) 


The y component in equation (3.87) may be rewritten as 

(3.90) 

Substituting equations (3.89) and (3.90) into equations (3.86), we may obtain r p . 

3.3.5 Principal Curvatures and Directions of the Pinion Cutter Surface at the Mean Contact 
Point 

The first principal direction of the pinion cutter surface at the mean contact point may be 
represented in the S p coordinate system as follows: 




sm r p = - 

sin A, 
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[ e Pj] - [- £ '/p(P)][- £ 'p(^m(^)][ i m(^ f; ( ^][ e p Jc ] 


(3.91) 


Using equation (2.9) and (3.91), we may obtain the first principal direction for the straight-edged 
cutter. It is 


e 


P U 


= ± 


sin sin r p 
cos T p 

cos S l sin r p 


(3.92) 


Using equations (2.17) and (3.91), we may obtain the first principal direction for the curve-edged 
cutter. The result is the same as for the straight-edged cutter, that is described in equation (3.92). 
In above equation, there are two senses. Only the direction which forms the smaller angle with 
the gear cutter first principal direction can be chosen. From the first order information we have 
already determined the parameter t p ; therefore, the first principal direction of the pinion cutter 
is also determined. The unit vector of the second principal direction of the pinion cutter surface 
may be obtained by rotating the unit vector of the first principal direction of the pinion cutter 
surface, e p ^ about the common normal, n, by an angle 7r/2. 

We use the concept discussed in Section A. 2 to derive the principal curvatures of the pinion 
cutter surface at the mean contact point. We recall that surfaces Ep and Ei are in line contact in 
the process of generation. Hence, using equation (A. 37), we obtain 

a n a 22~ a u= 0 ( 3 - 93 ) 


Substituting equations(A.31), (A. 32), and (A. 34) into (3.93), we obtain the first principal curvature 
of the pinion cutter 
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(3.94) 


K *I = 


K p n ( K i I cos2 a Pi + K i a sin2 C P i) - *!,*!„ 
• 2 2 
K B _ — K,_ sin <T PI — K, _ cos cr D1 


r i7 


u 


The second curvature of the pinion cutter is zero for a straight-edged cutter (see equation (2.12)) 
and =pl j R for a curve-edged cutter (see equation (2.20)). Since the principal curvatures and di- 
rections of the pinion cutter surface at the mean contact point have been determined, some data 
relating to pinion machine-tool settings may be obtained without any difficulty. 

Let us consider a straight-edged cutter first. Rewriting equation (2.10), we may obtain 


= 


1 

k P j tan xj) 


(3.95) 


We choose only the positive sign in equation (2.10) since we have specified the direction of the unit 
normal n. We may represent the mean contact point B in the S m (p) coordinate system as follows: 


[Bmi p )} = M m( p) p (P) M p(p)f [Bj] 


(3.96) 


where 



cos 6 X 0 sin 0 

0 1 0 T E m 

— sin 5, 0 costf, — L m 


0 0 0 1 


(3.97) 


and 
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-10 0 0 

0-100 

(3.98) 

0 0 10 

0 0 0 1 

Considering only the x component of the above equation, we obtain 

B rnx = -B fx cos S t + B h sin <5, (3.99) 

Using equation (2.38), we obtain r p as follows: 

r P = (B mx + u p cos^ p )tan^ P (3.100) 

Let us now consider the curve-edged cutter. Using equation (2.18), we obtain 




The cutter tip radius may be represented by 



(3.103) 
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3.4 Pinion Machine-Tool Settings 


There are five machine-tool settings rn pi , E m , 
the solution of this problem is the determination 
problem first. 


L rn , s p , and q p to be determined. The key to 
of the cutting ratio m pi . Let us consider this 


3.4.1 Determination of Pinion Cutting Ratio 

We will use the relations between principal curvatures and directions for the pinion cutter 
surface and the pinion surface to derive the pinion cutting ratio m pi . To apply the equations 
described in Section A. 2, we consider that surfaces Si and Sj r are equivalent, and that surfaces S p 
and S q are equivalent. Also, the following data are considered as given: (1) the principal curvatures 
of the pinion surface at the mean contact point, and K ln \ (2) the principal directions of the 
pinion surface at the mean contact point, e Xj and e ln \ (3) the coordinates of the mean contact 
point; (4) the unit normal at the mean contact point; (5) the coefficients a n , a 12 , and a 22 . 

The procedure to determine m F1 is as follows: 

^(ip) 

Step 1: Representation of u> 

The angular velocity of the pinion is represented by 

(3.104) 

The angular velocity of the pinion cutter is represented by 


-Ji) , (i) 

Ll> — ±UJ 


sin/i 1 

0 

I cos 
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cos 8 X 
0 

- sin 8 X 

Therefore, we may obtain the relative angular velocity u> as 

sin /ij — m pi cos 8 X 
0 

cos fi x + m pi sin 8 X 

Step 2: Representation of [ uj ne p ] 

The scalar [ £3 ne p J is represented by 

±(sin /ij - m pi cos 8 X ) 0 ±(cos ji x -f m px sin S x ) 

Tlj- TLy 72 ^ 

C Pj z 6p iy €P I x 

(3.107) 

— ^”1 [ ^ F; — n y e p, ) (^ x ^ p j y — ™y^Pi z ) ^ i j w pi 

+ [( n » e p Jt ~ n *Cp ) sin /x, + (n x e F ^ - n y e Pi J cos n, j } u> (1) 

= {c u m pi + c l3 )u> 

( 1 P) 

Step 3: Representation of [ u> ne p ^] 

( 1 P) 

The scalar [ uj ne p ^] is represented by 




_(p) , (i) 

u? — ±m pi u; 


(3.105) 


follows: 


(3.106) 
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-.(IP) ^ 
u> ne D 


(i> 


±(sin //j - m pi cos 6 X ) 0 ±(cos fi x + m pi sin 6 1 ) 




n 1t 


n 7 


= ± [( n y e p // - 


) S * n ( n ® € p jr/ ” n y € p n ) cos Ah W 


( 1 ) 


Step 4: Representation of V 


up) 


The velocity tr V { may be obtained by 


tr 


V 


( 1 ) 


£ (l> X B 


— ±U) 


(1) 


— By COS /Xj 

B x cos ji l - B z sin 
B y sin 


The velocity t 1 V [ ' may be obtained by 


tr 


V 


(P) 


AP) 


AP) 


u; x B T 0 fO Tn x u) 

(E m ± B y ) sin £, 

±(L m - B x sin £, - B y cos 6 l ) 
( E m ± By) cos 6 1 


(i) 

u) m 


The sliding velocity V (1P) is described by 


I (3.108) 
>) 


(3.109) 


(3.110) 
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Step 


fM lp ) 


\7 (p) 

= try “ fcrV 


= u> 


(i) 


ifB^cos^ - m pl (B m ± fly) Sin a J 

±(B X cos fi l — B z sin i n 1 ) m pi [L rn — (B x sin 6 t + B z cos 6 l )j 
±B y sin /ij - m pi (fl m ± By) cos b x 


( 3 . 111 ) 


. . ,r(lP] 1 T r( 1 F) 

5: Representation of v Pj and v Pjj 
Using equations (A. 33) and (3.107), we have 


a 


13 



(c n m 


PI + C 12J 


( 1 ) 


(3.112) 


Equations (A. 35) and (3.108) yield 


- C„U, 


(1) 


(3.113) 


From equation (A. 37) we have 


a n a 23 a i2 a i3 ® 


(3.114) 


Using equations (3.112) - (3.114), we obtain 


a i2 K P T Vpr ] a il K P rr ^ 


-(IP) _ 
P/7 " 


d" ( a n C 2 


^12 C 12 )] ^ 


( 1 ) 


(3.115) 


Moreover, we know that 
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(3.116) 


\ 


(ip) 



+ V}?? P 


II 


Considering only the x component in equations (3.111) and (3.116), we obtain 


V p) F) e Pli + Vpn e p Ul = [Tfiycos/i, - m pi (£ m ± B y ) sin tf, ] (3.117) 


Considering only the 2 component in equations (3.111) and (3.116), we receive 

Vp l F>e P h + Vp ! P e p I]; = [±B y sm^i l - m pi (E m ± B y ) cos 5, ] J 11 (3.118) 


Multiplying equation (3.117) by cos and equation (3.118) by sin S 1 , and adding the resulting 
equations, we obtain 


Vi 1P} 

v Pu 


T- 


By cos 7 t 


cos 6 j 


r IIz 


sin^, 


(i) , (i) 

-u; — t 4 uj 


(3.119) 


Substituting equation (3.119) into equation (3.117), we obtain 




Cn m I <lllKp n t * + a u c 22 a i2 C i2^ (i) , (l) 

~ m Pl + — W = (*!"»« +*»)“' 

K P, a i2 K Pl ) 


(3.120) 


Step 6: Representation of V * 1P * 

The matrix form of equation (3.116) may be represented by 
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(3.121) 


Step 


V' ,1P) 


+ KT‘„. 


V [IP] P 4-V ilp) e 

1 p i p K + V p n e p n v 


v?r<„. + ns’s, 


Iz 


Substituting equations (3.119) and (3.120) into equation (3.121), we have 


t(i p ) 


— uJ 


( 1 ) 


(t l m pi -f t 2 ) e Pj * ~^^4 c p nx 
(t x m p i -f t 2 )Cp ~ht 4 e p 

(*i m Pi + ^) e P /2 + ^p j/2 


(i) 


u n m pi + w ia 


U 21 ^P1 U 22 


L U 3i m P! + ^32 


( 3 . 122 ) 


7: Representation of [ nu;* lP) T r(l 5 ] 

The scalar [ nu; (1P) l r(1P) ] may be represented by 


^(1F) -(IF) 

nw V 






±(sin n l + m pi cos 6 X ) 0 ±(cos /i, + m pi sin 6 1 ) 


*11 Pi 1 12 


2 1 ""PI 1 22 


31 Pi 1 32 


/ 2 \ 

(i)l 

(f] m pi + V 2”V. + V Z ) 

u> 


(3.123) 


where 


v x = ± [(u n sin 6 X - it 31 cos 6 X )n y - ( n z cos 6, 4 - n x sin 6, )u 21 ] (3.124) 
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Step 


V 2 ~ =F [(^21 COS V 1 + U 22 S ^ n ^1 ) U X ~ (^21 S ^ n M 1 “ U 22 COS )™2 

- (i/jj cos /ij + u 12 sin^ + u 32 cos ^ - u 31 sin fi l )n y ) 


v 3 ~ T [( u 22 n x cos/Xj - (u 12 cos/q - u 32 sin /ijn y - u 22 n 2 sin /q] 


8: Representation of n • (u; (1) X * r V" <P) - u; (P} X * r F (1) ) 

— * ( pj 

The velocity t r V may be described by 


Substituting equations (3.109) and (3.122) into equation (3.127), we have 


T n(^> (i) 

t T \ — uJ 


— u ll m P1 T By cos /x, - tx 12 
— u 21 m Pl ^ B z sin /q ± B x cos /q — u 22 

~ U 3 i m Pi ± sin Mi - W32 


Vector (a; (I) 


x- tr F (P) ) 


is represented by 


-U) t t(p) 

a? X f r \ 



±[u 21 m P1 — ( B z sin /q — B x cos /q ± u 22 )] cos jq 
(T “u cos /x, ± tx 31 sin /x 3 )m pi - B y ^ u l2 cos fi 1 ± tx 32 
■f [ u 2i m pi - (B z sin /x, - £* cos /x, ± u 22 )] sin /x, 


(3.125) 

(3.126) 


(3.127) 

(3.128) 

Mi 

(3.129) 
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Vector (cJ (P) X t r V r(1) ) is represented by 


UJ 


x 


tr 


V 


t(l) 



— ( B z sin - B x cos p x )m pi sin 
— B y m pi sin 7^ 

-(5 z sin / u 1 - cos jj x )m pi cos 


Subtracting equation (3.130) from equation (3.129), we obtain 


-*(i) 

u; 


r 7( p > 


tr* 


P) 

UJ 


t(D 




^21 ^Pl ^22 


l ^3i d* h 32 


where 


/z n = ±u 21 cos /q — (I? x cos /ij — sin fi l ) sin 6 l 


h 12 = ( B z sm fi l — B x cos /q ± u 22 ) cos /q 


h 2J = By sin7 : u ll cos /x : ± u 31 sin/q 


h 22 = -{By ± u 12 cos/ij T u 32 sin/i,) 


/i 31 = u 3J S ^ n Pi — (-Bz: cos Hi — B z sin /ij ) cos 8 X 


h 32 = -(I? 2 sin /ij - 5 X cos /Xj ± n J2 )sin//, 


(3.130) 


(3.131) 


(3.132) 


(3.133) 


(3.134) 


(3.135) 


(3.136) 


(3.137) 
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Step 


Therefore, we may obtain n • (w* 1] X tr \ 


rAP) 


-*{P) 
— UJ 


x trV {1) ) as follows: 


n 




J i) 




— UJ 


J.P) 


trV W )={hm pl + /,) 


( 1)1 


<x> 


(3.138) 


where 


fi = n x h u + n y h 21 + n z h 31 (3.139) 


f 2 = n xh l2 + n y h 22 + n z h 32 (3.140) 


9: Representation of m Pi 

Using equations (A. 33), (3.107), and (3.120), the equation for a 13 may be represented by 


= - (k P j * 2 + C 1 2 ) 


( 1 ) 


UJ 


(3.141) 


Using equations (A. 35), (3.108), and (3.119), a 23 may be described by 





(3.142) 


Using equations (A. 36), (3.119), (3.120), (3.123), and (3.138), the expression for a 33 may be 
represented by 
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a 33 = [(2/Cp^,^ - v 2 - /j )m pi + + K Pjr t * - t> 3 - / 2 )] [a/* 1 ] 


(3.143) 


From equation (A. 39) we know that 


a i2 Q 33 fl 13 a 23 — ^ 


(3.144) 


Equations (3.141)-(3.144) yield 


l \2 {^Pj ^2 K P/J ^ 4 ^3 fl) { K Pjt 2 ^ C 12 )( K P J; ^4 C 22 ) 

a 12 (2K p ^tit 2 — v 2 — /j ) 


(3.145) 


3.4.2 Determination of parameters E m and Lm 


Parameters and Z m of the pinion machine-tool settings have been shown in Figures 12 
and 13. Since the pinion cutting ratio m pi has been determined, it is very easy to find these two 
parameters. We may determine vector V from equation (3.122). Applying equation (3.111), 
then, we obtain 


T B y cos fi 1 - v x 
m pi sin 6 X 


(ip) 




(3.146) 


L-rn — 


_ By cos - B z sin fi x V y 


(ip) 


m r 


+ B x sin 8 X + B z cos 8 X 


(3.147) 
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3.4.3 Determination of Pinion Radial Setting and Cradle Angle 


The determination of the pinion radial setting and the cradle angle is based on the consideration 
that the position vectors of the pinion tooth surface and head-cutter surface must coincide at the 
mean contact point. Equation (3.96) describes the mean contact point B in the S m {P) coordinate 
system. Considering the y and 2 components in equation (3.96), we obtain 


B 


( p ) 

y 


B j y T E-m 


(3.148) 


B (p) = Bj x sintfj + Bf z cos^ 1 - L m 

m z 


(3.149) 


For a straight-edged cutter, by using equations (2.38), (3.148), and (3.149), we have 


s p sin q p = ±Bj y + E rn ± u p sin^ sin r p (3.150) 


s p cos q p = B fx s'm6 1 + Bj : cos6 1 - Lm - UpS'mrppCOSTp (3.151) 


For a curved-edged cutter, by using equations (2.40), (3.148), and (3.149), we have 


s p sin q p = ±5y y + E m ± 


cos A p sin r p 


(3.152) 


SpCosqp = Bf t s'mS 1 + Bf z cos6 1 - Lm E 


cos A p cos T p 


(3.153) 
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Using sin q p -|- cos q p — 1, we eliminate q p and solve for pinion radius s p . Eliminating s py we 
may determine the pinion cradle angle q p . 
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CHAPTER 4 
CONCLUSION 

As it was mentioned in Chapter 1, the reduction of transmission errors of spiral bevel gears 
is a difficult problem. Although it is possible to generate conjugate spiral bevel gears, with zero 
transmission errors, we have to take into account that the gear are very sensitive to misalignment . 
Using the TCA programs we have found that even a small misalignment of gears results in discon* 
tinuity of functions of transmission errors that is accompanied with the jump of the function at the 
transfer points. Thus the idea of gears with non-zero transmission has to be complemented with 
the modification of the process for their generation that allows to reduce the sensitivity of gears to 
their misalignment. 

From the result of computation by TCA programs we know that gear misalignment causes 
a linear or almost linear function of transmission errors. Litvin has discovered that a sum of a 
parabolic function and a linear function represents again a parabolic function that is just translated 
with respect to the initial parabolic function. Then, if a parabolic function is predesigned, it 
becomes possible to keep the same level of transmission errors for aligned as well as misaligned 
gears. 

Gear misalignment is also accompanied with the shift of the bearing contact to the edge of gear 
tooth surface. To keep the shift of bearing contact in reasonable limits, it is necessary to limit the 
tolerances for gear misalignment and the respective value of predesigned parabolic function. 

In Chapter 2 the basic concept and methods of Gleason systems have been presented. Equations 
that describe the surface of the head cutter, which is either a cone surface or a surface of revolution, 
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have been derived. These equations covers the determination of position vectors, surface unit normal 

vectors, principal curvatures, and principal directions. 

Mathematical models for geometry of spiral bevel gears have been also proposed in Chapter 2. 
The gear surface is represented as an envelope of the family of the tool surfaces. The tool surface 
and being generated gear surface are considered conjugate ones. Based on the geometric properties 
of conjugate surfaces, the equation of meshing has been derived. 

The determination of pinion machine-tool settings is based on the method of local synthesis. 
The first derivative of gear ratio, the tangent to the contact path, and the dimensions of the contact 
ellipse of the gear surface at the mean contact point are considered as input to local synthesis. Thus 
the level of transmission errors and the bearing contact are under control. It provides the optimal 
conditions of meshing for the gear surfaces being in meshing at, and within the neighborhood of, 

the mean contact point. 

Equations that determine the principal curvatures and directions at the mean contact point on 
the pinion surface have been derived. They are functions of the principal curvatures and directions 
at the mean contact point on the gear surface and the input of local synthesis. Based on the 
information on the characteristics of the pinion surface of the zero (position), first (normal), and 
second (principal curvatures and directions) orders, equations that determine the pinion basic 
machine-tool settings have been derived. 

In Appendix A the basic concept and methods of theory of gearing that have been used in this 
work have been presented. Numerical examples are given in Appendix B. These examples include 
determination of machine-tool settings and results of computation by TCA programs. Computer 
programs have been developed, that include machine-tool settings and TCA. They are listed in 
Appendix C. The computer programs cover determination of machine-tool settings for straight - 
lined as well as curved blades. The developed TCA programs allow to simulate the meshing of 
aligned and misaligned gears. 
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APPENDIX A 

GEOMETRY AND KINEMATICS OF GEARS IN THREE DIMENSIONS 


A.l Concept of Surfaces 1 

Most of the ideas underlying gear theory are based on strict definitions proposed in the field 
of differential geometry. In what follows we introduce the concept that is applied in this report. 

All in all we require that our functions can be differentiated at least once and usually more times. 
Accordingly we say a function F belongs to class C n on an interval 1 if the nth order derivative 
°f F exists and 1S conti nuous on 1. In addition, we denote the class of continuous functions by C°. 

A parametric representation of a surface E is a continuous mapping of an open rectangle 3?, 
given in the plane P of the parameters (u,t>), onto a three-dimensional space £ 3 such that 

B(u,v)€C°, (u,v)e3? (A.l) 

where B is the position vector which determines the point surface (Figure 21). The vector func- 
tion B(u,v) may be represented by 

B(u,v) = B x (u,v)T+ B y (u,v)j+ B z (u,v)k (A. 2) 


'Adopted from the manuscript of the book “Theory of Gearing” by Litvin, in press by NASA. 
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where z, j , and k are unit vectors of the coordinate axes. 

We call a surface point B(u , v) a regular point if at this point 


B u x B v ^ 0 


(A.3) 


where 


B u 


dB 



d£ 

dv 


A surface is called a regular one if each point on it is a regular point. 

A regular surface has the following properties: 

• It is at least class of C 1 . 

• There is a one-to-one correspondence between the points of plane P (of the parameters (u,v)) 
and the three-dimensioned space B 3 . 

• A regular surface has a tangent plane at all its points. 

The normal vector N to the surface at a point B is 


N = B u x B v 


(A- 4) 


and its unit normal is represented by 


N _ Bu X By 

|A| “ B u x B v 


(A.5) 
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The direction of the surface normal N and unit normal n, with respect to the surface, depends on 
the order of the factors of the cross product (equation (A. 4)). By changing the order of the factors, 
we may change the direction of the normal to the opposite direction. 

A surface is uniquely determined by certain local invariant quantities called the first and second 
fundamental forms. The first fundamental form of a surface is defined by 

I — dB ■ dB = (B n du + B v dv) ■ (B u du + B v dv ) 

= (B u • B u )du + 2 (B u • B v )du dv + (B v • B v )dv (A.6) 

= E du + 2 F du dv + G dv 2 


where we set 


E — B u ■ B u , F = B u - B v , G — B v B v 


The second fundamental form is 

II — - dB • dn = ~{B U du + B v dv) • (n u du + n v dv) 

= -{B u ' n u )du - (5 U • n v + B v • n u )dudv - (B v ■ n v )dv (A. 7) 

— L du -f 2 M du dv + A 7 di; 2 


where we have 


L — B u • , 


AI — { B u ’ ti v B v • N — B v • ti v 


The second fundamental form exists only if the surface is at least class C 2 . In this report we will 
consider all the gear tooth surfaces as regular surfaces with class at least C 2 . 

c- 2- 
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On a given surface various curves pass through a common point B and have the same unit 
tangent vector f at B (Figure 22). One of these curves (designated by L 0 ) is located on the 
plane P, which is drawn through the unit tangent vector f and the surface unit normal n. The 
curvature of curve L 0 is called normal curvature. Since the unit tangent vector f of the surface 
may have different directions on the surface, for each direction there is a normal curvature. The 
normal curvature is a function of the first and second fundamental forms: 

_ II _ Ldu +2M dudv + N dv 

Kn ~ J ~ £du 2 +2Fdtxdu + Gdu 2 ( A ' 8 ) 


The extreme value of the normal curvature taken at a certain point of the surface are called the 
principal curvatures. The directions of the normal sections of the surface with the extreme normal 
curvatures are called the principal directions. Equation (A. 8) yields 


IF — n n (E du + 2 F du dv + G dv ) — (L du + 2 M du dv + N di ? ) = 0 (A. 9) 


For a given point on the surface, E, F, G, L, M, and N are constant. The normal curvature K n 
is a function of the ratio du and dv. Therefore, equation (A. 9) is an identity of du and dv. From 
calculus, the partial derivative 


dF 

d du 


= 0 


(A. 10) 


Substituting equation (A. 9) into equation (A. 10), it yields 

n n (Edu+ F dv) - ( Ldu + M dv) + -^~{E du + 2 Fdudv + G dv 2 ) = 0 (A. 11) 
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Also, the partial derivative 


d_T_ 

d dv 


= 0 


(A. 12) 


Substituting equation (A. 9) into equation (A. 12), it yields 


K n{E du -+■ G dv ) — (M du + N dv ) T 


9Kji 
d dv 


E du + 2 F dudv + G dv = 0 


(A. 13) 


Recall that the principal curvatures are the extreme values of the normal curvature K n . Thus 
dn Tl /ddu = 0 and dn n jddv = 0 if K n is the principal curvature. Equations (A. 11) and (A. 13) yield 


(KnE — Z) du + (/c n F — Af) dv — 0 


(A. 14) 


and 


{n n F - M) du + ( k h G - N) dv = 0, 


(A.15) 


respectively. Solving the homogeneous system of equation (A.14) and (A.15) by eliminating du 
and dv , we obtain 


{EG - F 2 ) K 2 n - (EN - 2 FM + GL)K n + (LN -M*) = 0 


(A.16) 


The discriminant of equation (A.16) is 
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(A. 17) 


A = (EN -2FM + GL) 2 - \(EG - F 2 ){LN - M 2 ) 


2 

= \(EN - GL) - 2£{EM - FL)\ + *( EG ~ F --^( EM - FL ) 2 


Equation (A. 17) shows that the discriminant is greater than or equal to zero. Thus the equation 
has either two distinct real roots— the principal curvatures at a nonumbilical point, or a single real 
root with multiplicity two — the curvature at an umbilical point. The discriminant is equal to zero 
if and only if 


EN - GL = EM - FL = 0 


(A.18) 


Since E / 0 and G / 0, equation (A.18) can be shown to be identically equal to 

L ----- n 

E~ ~F ~ G ~ 


(A.19) 


Considering equations (A. 8) and (A.19), we obtain 

K n — M 


(A. 20) 


This means that the principal curvature is the same as the normal curvature at any direction. 
Thus each direction may be considered as a principal direction. Any point which is on a plane or 
at which a surface turns into a plane 2 is an umbilical point. Any point on a spherical surface 3 is 
also an umbilical point. 

Two distinct principal curvatures can always be obtained at a nonumbilical point. These two 

curvatures correspond to two distinct principal directions. By canceling /c n from equations (A. 14) 
2 The normal curvature on each direction is zero. 

3 The normal curvature on each direction is the inverse of the radius. 
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and (A. 15), we have the following equation for principal directions 


{EM - FL ) du - ( GL - EN) dudv + (FN - GM) dv = 0 (A. 21) 


The discriminant of the above equation is identical to equation (A. 17). At a nonumbilical point 
equation (A. 21) can be represented as a product of two co-factors (A z du + B z dv)(i = 1,2) since 
the discriminant is larger than zero. This means that it represents two perpendicular directions. 
Thus we may conclude that at a nonumbilical point there exist two distinct principal curvatures in 
two perpendicular directions. 

After representing of E, F, G, L, M , and A r in the form of B u , B v , n u , and n v by equations (A. 6) 
and (A. 7), equations (A. 14) and (A. 15) will yield 

B u ■ (K n dB + dn) = 0 (A. 22) 

B v ■ (k ti dB + dn) = 0 (A. 23) 


Obviously, 


n ■ (k ti dB + dn) - 0 


(A. 24) 


Therefore, K n dB + dn is a zero vector since it is orthogonal to B u , B v , and n. In short, we have 


dn = 


— K n dB 


(A. 25) 
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The above equation, which completely characterizes the principal curvatures and directions, is 
called Rodrigues' formula. This formula simplifies the calculations to obtain principal curvatures 
and principal directions. The matrix form of Rodrigues' formula is 


dn x , dn x 

--du+^-dv 

du dv 

3 Tty . dTiy 

— V. du + -r-2- dv 
du dv 

dn z dn z 

— du + dv 
du dv 


~ K t.n 


dB x 

du 

dBy 

du 

dB z 

du 


du + 
du + 
du + 


dB x 

dv 

dBy 

dv 

dB z 

dv 


dv 

dv 

dv 


(A. 26) 


Matrix equation (A.26) yields three scalar equations in three unknowns, the ratio du/iv, and the 
principal curvatures «, and K„. Using any two of the scalar equations we may develop a quadratic 

equation (provided dv ± 0) 


A 2 



+ Al p + A D = 0 

dv 


(A.27) 


The two roots of this equation correspond to two principal directions on the surface. By putting 
both roots into the third scalar equation, we may determine the principal curvatures k, and 

It is possible to have either positive or negative principal curvatures. The sense of the principal 
curvature depends on the location of the center of curvature on the normal. The principal curvature 

is positive if the center of curvature is located on the positive normal. 

The normal curvature on each direction may be expressed in terms of principal curvatures. This 

is so called Euler’s Theorem. That states 

Kn = KjCos & + & (A. 28) 

where w is the angle formed by the tangent to the normal curvature and principal direction with 
curvature K r 
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A ' 2 j*;gI atlons Between Principal Curvatures and Directions for Mating Surfaces 


Consider two gear surfaces Ejr and Eg which are in meshing. Moreover, we have the following 
assumptions: 

1. The rotation angles, <f> T and <p Q , of both gears are given; 

2. The function (j> Q {(f)^) has continuous derivatives of second order; 

3. The angular velocity of gear T is constant. 

Then relations between principal curvatures and directions of these mating surfaces may be deter- 
mined. Such relations were first proposed by Litvin [12] and then extended for the case m' ^ 0 
by Litvin and Gutman [3], where m TQ = / J a) i s the gear ratio. 

The relations may be expressed by a system of three linear equations in two unknowns r V^ 
and r V%: 

a 3 i rV, a/ + a i2 t VqJ = a j3 ( j = 1,2,3) (A. 29) 


where r V 0/ and are the projections of the relative velocity r V tn at the contact point B on 

the principal directions on surface Eg. The equation may be represented by a symmetric augmented 
matrix [A]. That is 


a il a i2 a i3 


[A] 


a 21 a 22 fl 23 


a 31 fl 32 fl 33 


(A. 30) 


Here 
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= K _ - K _ COS <7 

C/ *7 


k T sin <r — 
T u 




k t - k t 

^ " - yj cos 2 <t (A. 31) 


= = 


k Tj ~ K r 


— sin 2a 


(A. 32) 


— k 0 K - kj ne n 

C 7 Qj l «/ J 


(A. 33) 


a 00 — K _ — K_ Sin <7 — K_ COS <7 = 

22 Ci7 ^7 


K _ -f 
_^_i T n 


k t, - K r n 


Q-u 


cos 2(7 (A. 34) 


a 


23 


— K a 
£j7 


y[FQ) 
<2 77 


r -(^C) _ 

[a> ne 


C/7 


(A. 35) 


a 


33 



nu; V 


(A. 36) 


n 


U> 


X trV 


(C) 


ot 

— u; 


trV 


i7<* 


) - (O' 


m , 


Q.T 


(n 


k Q ) • (B - 0 


*0 T ) 


k t ^ and k 7 ^ are the principal curvatures at the contact point B of gear T , 

K Qr and K Qn are the principal curvatures at the contact point B of gear Q, 

e Qi and e Q ^ are the unit vectors of the principal directions at the contact point B of gear Q, 
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a is the angle measured counterclockwise from e T , the unit vector of the principal direction at 
the contact point B of gear T , to e Q , 
c? (r) and u} {Q) are the angular velocities of gears T and Q, respectively, 
uJ (;rC) is the relative angular velocity of gear T with respective to gear Q, 
n is the common unit normal vector, 

V {TQ) is the relative velocity of the contact point on gear T with respect to the same contact 
point on gear Q, 

V*/* and v£® are the projections of V 1 *® on the e Qj and e Qj] , respectively, 

and fr F (C) are the transfer velocities of contact point B on gear T and gear Q, respectively, 
B is the position vector of the common contact point B : 

0 Q 0 T is the position vector from 0 Q to 
k is the unit vector of the axis of rotation of gear Q, and 

m! is the derivative of the rotation ratio of gear Q to gear T . It is represented as 


m 


/ 

QF 


d 

d<t> T 


rn 


QT 




where <\> T is the rotation angle of gear T , and 


(Q) 

w 

_ (?) 
w 


Totally, there are two cases of tangency of gear tooth surfaces: 

1. The surfaces and are in line contact and B is just a point of the instantaneous line of 
contact. 
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2. The surfaces E T and H Q are in point contact and B is the single point of tangency at the 
considered instant. 

In the case of line contact of mating surfaces, the rank of matrix [^4] is equal to one. Thus all 
determinants of the second order formed from the elements of [.4] are zero. This yields 



_ G 12 

— Ql3 

^13 

1 

kP | 

a 33 


a i2 

_ ^22 

— ^23 

<*13 

°23 

a 33 


Using equations (A.31)-(A.39) we obtain 


tan 2cr 


- a 


^ Q 13 fl 23 
- ( k Q/ - 


2 11 


K 


(A.37) 


(A. 38) 


(A. 39) 


(A. 40) 


K r — — 

T i T n 


- Q 13 +( K gj 

a 33 cos 2cr 


(A.41) 


2 2 
a , *f 


r ”f‘ K-t — ( K Q -f )— “ 

r ; T n ' S/ Qn' a 


(A.42) 


For the case when surfaces Ejr and Eq are in point contact, the rank of matrix [A] described 
by equation (A. 30) is two. Consequently, 
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= 0 


Equation (A. 43) yields the relation 


5 K T n ’ k Qj i K Q n ’ a ) ~ ^ 


(A.43) 


(A. 44) 


In general, equations of the generated surface are evidently much more complicated than those 
of the generating one. Therefore, a direct way to obtain the principal curvatures and directions of 
the generated surface is a very difficult task. This work can be significantly simplified if we apply 
the relations, described in this section, between principal curvatures and directions of meshing 
surfaces. 


A. 3 Relative Normal Curvature 

The relative normal curvature, K r , of two mating surface, Ejr and Eg, at the contact point B 
is defined as the difference of the normal curvatures of both surfaces taken in a common normal 
section of surfaces and represented as 


(C) (;r) / A 
K r = Kn - *„ (A. 45) 

Suppose the common normal section form an angle w with the unit vector e Q and an angle (zj + 
<j) with the unit vector e T (Figure 23). According to Euier’s Theorem (equation (A. 28)), we obtain 

(Q) 2 2 (;F) 2 2 

K n =k Qi cos & + k Qij sin w K n - k Tj cos (w + a) + sin (ttr + cr) (A. 46 ) 
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Substituting equation (A.46) into (A.45), after simple transformations we get 


K r = (k Qi - cos 2 a - Kj: n sin 2 a) cos 2 w + (n Q[j - sin a - cos a) sm lv 

+ -(it- - k t ) sin2crsin2tx7 
' o ' T i T u } 


(A. 47) 


Equation (A. 47) and expressions for a n , a 12 , and a 22 in equations (A. 31), (A. 32), and (A. 34) yield 


k t = l - [a n + a 22 + (a u - a 22 )cos2CT] + a 12 sin2^ 


(A. 48) 


The extreme values of function k t (zzt) may be determined by 




(/c r ) = 0 


(A. 49) 


Thus we obtain 


2a 2 

tan 2 w — 


(A. 50) 


This equation has two solutions w, and zu 2 . Moreover, \w 1 - w 2 \ - 7r/2. This means that there are 
two perpendicular directions for the extreme relative normal curvatures. Using equations (A. 48) 
and (A. 50) the extreme values of the relative normal curvatures are represented by 


Kr ~ 2 


( a n + a 2 2 ) ± \/( a 11 - a 2 2 ) 2 + 4 a i2 


(A.51) 


We may determine whether or not two surfaces interfere each other by the concept of relative 
normal curvatures. If two surfaces contact at a point with any interference, the sign of the relative 
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normal curvature in each direction must remain the same. In other words, the product of two 
extreme values of the relative normal curvatures is positive. Equations discussed in this section 
were first proposed by Litvin [9]. 


A. 4 Contact Ellipse 

In theory the tooth surfaces of a pair of spiral bevel gears are in contact at a single point at 
every instant. In practice the surface of the solids is deformed elastically over a region surrounding 
the initial point of contact, thereby bring the two bodies into contact over a small area in the 
neighborhood of the initial contact point [13, 14]. Such an area is an ellipse whose center of 
symmetry is the theoretical point of contact and the dimensions depend on the elastic approach 
and principal curvatures and directions of the contacting surfaces. If the approach of surfaces under 
the action of load is given, the size and orientation of the contact ellipse can be defined as a result 
of a geometric solution. Litvin[9, 15] has investigated the mathematical modeling of the contact 
ellipse. 

Let us now consider that two surfaces and £2 are contact at a single point B . The principal 
curvatures, and k x of £1 and k 2j and k 2 of £ 2 , at point B are known. Also known are unit 
vectors e Xj and e l ^ , which are directed along the principal directions of £1 at point B , and e 2 
and e 2jj , which are directed along the principal directions of £ 2 at point B . Unit vectors e lj and e 2 ^ 
determine the tangent plane (Figure 24). Angle <r 12 , which is measured counterclockwise from e x 
to e 2 ^ is also determined since e lj and have already been known. Then the contact ellipse may 
be described as 



(A. 52) 


in which £ and 77 are coordinates with respect to the £ and 77 axes with origin at the contact point B . 
The lengths of semiaxes a and b are 
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c 



(A. 53) 



b = 




£ 

B 


where £ is the approach, and 




2k 1a k 2A cos 2<t 12 



(A. 54) 




^2A cos 2 <t 12 



(A. 55) 


where 


= «1, +*1„- 


- k 2j + k 2 


(A. 56) 


K 1A “ K l; ^ 1 jj 1 


K A J-k 

2 I 2// 


(A. 57) 


The angle a, which determines the orientation of the ellipse may be obtained by equations 


cos 20 j 


k 2A cos 2<t 12 


- 2k iaSa cos2ct 12 + k 2A 


(A. 58) 


and 


sin 2a x 



2 k, 2A sin 2cr 12 
2k ia K 2A COs2(7 12 + K 2A 


(A. 59) 
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Finally 


(A. 60) 

Note that the angle a l is measured counterclockwise from the rj axis to the unit vector e 
Since 

the semimajor axis of the contact ellipse may be determined by the following conditions: 

• The length of the semimajor axis, which is along the rj axis, is b 
if n 2Tt > k 1s or \A\ > \B\. 

• The length of the semimajor axis, which is along the ( axis, is a 
if k 12 > k 2S or \B\ > \A\. 



arctan 


sin 2 Qj 
1 + cos 2a 
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APPENDIX B 

NUMERICAL EXAMPLES 

In this section, we will use the synthesis method discussed in Chapter 3 to determine the 
machine-tool settings for a pair of spiral bevel gear drive, and then we will use the TCA to simulate 
the meshing of this pair under alignment and misalignment conditions. Two cases are considered 
here. Both cases use straight blades to cut gears, but for the pinion, case 1 uses straight blades, 
and case 2 uses curved blades. 

The major blank data is represented in Table 2. Table 3 shows the input for case 1, and Table 4 
shows the input for case 2. The output for the gear machine-tool settings is shown in Table 5, 
which is the same for both cases. For the pinion machine-tool settings, case 1 is shown in Table 6, 
and case 2 is shown in Table 7. 

Two conditions of misalignment are considered when the TCA is applied to simulate the mesh- 
ing. They are the shift of pinion along its axis, which is denoted by A A, and the error of pinion 
shaft offset, which denoted by AV. We consider that A A is positive when the mounting distance 
of pinion is increased. The sense of AV is the same as y } shown in Figure 18. The output of the 
TCA is shown from Figure 25 to Figure 34 for case 1, and from Figure 35 to Figure 44 for case 2, 
respectively. 
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Table 2: BLANK DATA. 



Pinion 

Gear 

Number of Teeth 

10 

41 

Diametral Pitch 

5.559 

Shaft Angle 

90° 

Mean Cone Distance 

3.226 

Outer Cone Distance 

3.796 

Whole Depth 

0.335 

Working Depth 

0.302 

Clearance 

0.033 

Face Width 

1.139 

Root Cone Angle 

12°T 

72°25' 

Mean Spiral Angle 

— 

35° 

Hand of Spiral 

R.H. 

L.H. 


Table 3: INPUT DATA FOR CASE 1. 



Gear Convex Side 

Gear Concave Side 

Gear Blade Angle 

20° 

Gear Cutter Average Diameter 

6 

Gear Cutter Point Width 

0.08 

First Derivative of Gear Ratio 

-0.0035 

0.0052 

Semimajor Axis of Contact Ellipse 

0.171 

0.181 

Contact Path Direction Angle 

90° 

75° ~ ’ 
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Table 4: INPUT DATA FOR CASE 2. 


Gear Convex Side 

Gear Concave Side 

Gear Blade Angle 

20° ! 

Gear Cutter Average Diameter 

6 i 

Gear Cutter Point Width 

0.08 

First Derivative of Gear Ratio 

-0.0037 

0.0055 

Semimajor Axis of Contact Ellipse 

0.171 

0.171 

Contact Path Direction Angle 

90° 

cn 

o 

Radius of Blade 

40. 

50. 


Table 5: GEAR MACHINE-TOOL SETTINGS. 


Radial 

2.87798 

Cradle Angle 

58.6365 

Ratio of Roll 

0.973748 




Table 6: PINION MACHINE-TOOL SETTINGS WITH STRAIGHT BLADE. 



Pinion Concave Side 

Pinion Convex Side 

Blade Angle 

16.5561° 

22.9907° 

Tip Radius of Cutter 

2.96469 

3.07037 

Radial 

2.99331 

2.69783 

Cradle Angle 

63.1869° 

54.1910° 

Ratio of Roll 

0.22900 

0.25348 

Machining Offset 

0.17404 

-0.24459 

Machine Center to Back + 
Sliding Base 

0.021231 

0.052118 


Table 7: PINION MACHINE-TOOL SETTINGS WITH CURVED BLADE. 



Pinion Concave Side 

Pinion Convex Side 

Blade Angle 

16.5561° 

22.9907° 

Blade Center 

(11.557, 0., -35.309) 

(19.685, 0., 49.006) 

Tip Radius of Cutter 

2.98467 

3.04386 

Radial 

2.95578 

2.74261 

Cradle Angle 

63.0025° 

54.0900° 

Ratio of Roll 

0.23157 

0.24915 

Machining Offset 

0.12042 

-0.18825 

Machine Center to Back + 
Sliding Base 

0.01690 

0.03605 
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GEflR CONVEX SIOE 
AUSwCNT 




Figure 25: Straight-edged blade, gear convex side, alignment. 
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GEAR CONCAVE SIOE 
flUGWCNT 




Figure 40: Curved-edged blade, gear concave side, alignment. 
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appendix c 


LISTING OF COMPUTER PROGRAMS 


* it * * * * A Vr * * * * * * s'c * * * * * * * * 3'c * >V * * * * * * * 3V * * * * 3V * * * * Vr * * * * 3V * * * * * 3 *c * 3V 3V * Vr 3’c 3V * 3V * * * 3 V * sc 3 f 

* 


* Gleason's Spiral Bevel Gears 

* * 

* Basic Machine-Tool Settings and Tooth Contact Analysis 

* * 

* Straight Blade to Cut the Pinion * 

3’C 

3? 

* * * * * * * * * * * * * * * * * * * * * * * * * * j't * * A * * * * * * * * Vc * * * * Vc * * * * * * * * * jV sV * * * * * * * * X * 5f 5C * w * * 


IMPLICIT REAL*8(A-H,K,M-Z) 

REAL *8 X (1) , F (1) ,FI(1) , PAR (6) ,LM,TX(5) ,TF(5) ,TF1(5) ,TPAR(19) , 
AZSP (1,1) , WORKP (1) ,AZS(5,5) ,WORK(5) 

CHARACTER *8 HG, HNGR 
DIMENSION IPVTP (1) , IPVT (5) 

EXTERNAL PCN,TCN 
COMMON /P 1 /PAR 
COMMON / T 1 / TP AR 
COMMON /AO /HG 

COMMON/Al/pll ,pl2,pl3,p21 ,p22,p23,p31 ,p32,p33,pl,p2,p3 

COMMON/A2/SG , CSRT2 , QG, SNPSIG , CSPSIG 

COMMON/ A3 /TND1 , TND2,RITAG 

COMMON/A4/CSD2 , SND2 , CSPIT2 , SNPIT2 

COMMON / A5 / C S QG , SNQG, THETAG 

COMMON/Bl/CSPHl 1 , SNPHl 1 , SP , EM, LM, CSRT1 ,CSDl , SNDl , CSPSIP , SNPSIP 

COMMON/B2/CSPIT1 , SNPIT1 ,MPl ,MG2 ,QP 

COMMON/B3/B2fx,B2fy,B2fz 

COMMON/B4/CSPH2 , SNPH2 , CSPH2 1 , SNPH2 1 

COMMON/C 1 /UG , CSTAUG , SNTAUG 

COMMON/C2/N2fx,N2fy,N2fz 

COMMON/DI /UP , CSTAUP , SNTAUP 

COMMON / E 1 / XB f , YBf ,ZBf 

COMMON/ Fl/PHIGO 

COMMON/Gl/DAl ,DVl 


* INPUT THE DESIGN DATA 


3V 

* TNI : number of pinion teeth 

* sec. 3.1 
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* TN2 : number of gear teeth 

* sec, 3.1 

* RTldg, RTlmin : root angle of pinion (degree and arc minute, respec - 

* tively) 

* sec. 3.1 

* RT2dg, RT2min : root angle of gear (degree and arc minute, respec- 

* tively) 

* sec. 3.1 

* SHAFdg : shaft angle (degree) 

* sec. 3.1 

* BETAdg : mean spiral angle (degree) 

* sec. 3.1 

* ADIA : average gear cutter diameter 

* sec . 3.1 

* W : point width of gear cutter 

* sec. 3.1 

* A : mean cone distance 

* sec. 3.1 

* ALPHdg : blade angle of gear cutter (degree) 

* sec. 3.1 

* DLTXdg : angle measured counterclockwise from root of gear to 

* the tangent of the contact path (degree) 

* gear convex side 

* fig. 19 

* DLTVdg : angle measured counterclockwise from root of gear to 

* the tangent of the contact path (degree) 

* gear concave side 

* fig. 19 

* M21XPR : first derivative of gear ratio 

* gear convex side 

* sec , 3.1.1 

* M21VPR : first derivative of gear ratio 

* gear concave side 

* sec. 3.1.1 

* AXILX : semimajor axis of contact ellipse 

* gear convex side 

* eq. (3.76) 

* AXILV : semimajor axis of contact ellipse 

* gear concave side 

* e q. (3.7 6 ) 

* HNGR : hand of gear ( 1 L 1 or 1 R * ) 

* DA : amount of shift along pinion axis 

* + : pinion mounting distance being increased 

“ : pinion mounting distance being decreased 

* DV : amount of pinion shaft offset 

* the same sense as yf shown in fig. 18 

* DEF : elastic approach 

* eq. (3.76) 

* EPS : amount to control calculation accuracy 

* 

* OUTPUT OF THE BASIC MACHINE-TOOL SETTINGS 

* 

* PSIGdg : gear blade angle 
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* PSIPdg 

* RG 

* RP 

* SG 

* SP 

* QGdg 

* QPdg 

* MG2 

* MPl 

* EM 

* LM 


pinion blade angle 

tip radius of gear cutter 

tip radius of pinion cutter 

gear radial 

pinion radial 

gear cradle angle 

pinion cradle angle 

gear cutting ratio 

pinion cutting ratio 

machining offset 

machine center to back + sliding base 


DATA TNI , TN2/ 10. D00, 41 . D00/ 

DATA RTldg,RTlmin/12.D00, 1 . D00/ 
DATA RT2dg , RT2min/72 . D00, 25 ,D00/ 
DATA SHAFdg , BETAdg/90 . D00 , 35 . D00/ 
DATA ADIA/6.ODO0/ 

DATA W/0.08D00/ 

DATA A/3.226D00/ 

DATA ALPHdg/20 . D00/ 

DATA DLTXdg/ 90.D00/ 

DATA DLTVdg/ 75.D00/ 

DATA M21XPR/-3.5D-03/ 

DATA M21VPR/5.2D-03/ 

DATA AXILX/0. 1710D00/ 

DATA AXILV/0. 1810D00/ 

DATA HNGR/'L'/ 

DATA DV , DA/ 0 . D00 , 0 . D00/ 

DATA DEF/0. 00025D00/ 

DATA EPS/l.D-12/ 


DAl=DA 
DV 1=DV 
HG=HNGR 

* 

* CONVERT DEGREES TO RADIANS 

* 

CNST=4 . D00*DATAN ( 1 . D00) / 180 . D00 

RITAG=90 . D00*CNST 

DLTX=DLTXdg*CNST 

DLTV=DLTVdg*CNST 

RTl=(RTldg+RTlmin/60.D00) *CNST 

RT2= (RT2dg+RT2min/60 . D00) *CNST 

BETA=BETAdg*CNST 

PSIG=ALPHdg*CNST 

SHAFT=SHAFdg*CNST 

CSRT2=DCOS (RT2) 

SNRT2=DSIN(RT2) 

CSRTl=DCOS (RT1) 

SNRT1=DSIN (RT 1 ) 
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* CALCULATE PITCH ANGLES 

Vc 

MM21=TNl/TN2 
c eq. (3.1) 

PITCH2=DATAN (DSIN (SHAFT) / (MM21+DC0S (SHAFT) ) ) 

IF (PITCH2 .LT. O.DOO) THEN 
PITCH2=PITCH2+18O.DO0 
END IF 

CSPIT2=DCOS (PITCH2) 

SNPIT2=DSIN (PITCH2) 

c eq. (3.2) 

PITCH1=SHAFT-PITCH2 
CSPIT1=DC0S (PITCH1) 

SNPIT1=DSIN (PITCH1) 

* 

* CALCULATE DEDENDUM ANGLES 

* 

c eq. (3.3) 

D1=PITCH1-RT1 

D2=PITCH2-RT2 

CSD1=DC0S(D1) 

SND1=DSIN(D1) 

TND1=SND1/CSD1 
CSD2=DCOS (D2) 

SND2=DSIN(D2) 

TND2=SND2/CSD2 

* 

* CALCULATE GEAR CUTTING RATIO 

* 

c eq . (3.7) 

MG2=DSIN (PITCH2) /CSD2 

* 

* FOR GEAR CONVEX SIDE 1=1, FOR GEAR CONCAVE SIDE 1=2. 

■>v 

DO 99999 1=1,2 
IF (I .EQ. 1)THEN 
WRITE (72,*) 'GEAR CONVEX SIDE' 

DLTA=DLTX 
M2 1 PRM=M2 1 XPR 
AXIL=AXILX 
ELSE 

WRITE (72,*) 'GEAR CONCAVE SIDE' 

DLTA=DLTV 
M21PRM=M21VPR 
AXIL'AXILV 
END IF 
WRITE (72 , *) 

c eq. (3.76) 

AXIA=DEF/ (AXIL*AXIL) 

* 

* CALCULATE GEAR BLADE ANGLE 

* 

c sec. 2.2 
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IF (I .EQ. 2) THEN 
PSIG=180.D00*CNST-PSIG 
END IF 

CSPSIG=DCOS (PSIG) 
SNPSIG=DSIN (PSIG) 
CTPSIG=CSPSIG/SNPSIG 

* 

* CALCULATE CUTTER TIP RADIUS 

* 

c eq . (3.8) 

IF (I .EQ. 1)THEN 
RG= (ADIA-W) / 2 . DOO 
ELSE 

RG= (ADIA+W) / 2 . DOO 
END IF 


* CALCULATE RADIAL 

Vc 


— eq. (3.9) 

IF (I .EQ. DTHEN L 

SG=DSQRT (ADIA ,,c ADIA/4.D00+A*A*CSD2 Vf CSD2-A x ADIA ,c CSD2 x DSIN (BETA) ) 


* CALCULATE CRADLE ANGLE 

* 

C QG=DACOS( (A*A*CSD2*CSD2+SG*SG-ADIA*ADIA/4.D00) / (2 . D00*A*SG*CSD2) ) 

CSQG=DCOS (QG) 

SNQG=DSIN (QG) 

END IF 

A 

PAR (1) =RG*CTPSIG*CSPSIG 
PAR(4)=RG*CTPSIG 


* CALCULATE PHIG 

PHIG=O.DOO 
PHIGO=PHIG 
CSPHIG=DCOS (PHIG) 

SNPHIG=DSIN (PHIG) 

* 

IF(HG .EQ. 1 L 1 ) THEN 
IF (l .EQ. 1)THEN 

* Mmc=Mms*Msc 

c eq. (2.26) 

CALL COMB I (ml 1 , ml 2 , ml 3 , m2 1 , m22 , m23 , m3 1 , m32 , m33 , ml , m2 , m3 , 

1 .DOO, 0 .DOO, O .DOO , 0 .DOO , CSPHIG, SNPHIG, O.DOO, -SNPHIG, CSPHIG, 

0. DOO, 0. DOO, 0. DOO, 

1 . DOO , 0 . DOO , 0 . DOO , 0 . DOO , CSQG , -SNQG , 0 . DOO , SNQG , CSQG , 

0 . DOO , -SG*SNQG , SG*CSQG) 

END IF 

* Mpc^Mpm^Mmc 

c eqs. (2.25), (3.13) 

CALL COMBI (pll ,pl2,pl3,p21 ,p22,p23,p31,p32,p33,pl ,p2,p3, 
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. CSD2 , 0 . DOO , -SND2 , 0 . DOO , 1 . DOO , 0 . DOO , SND2 , 0 . DOO , CSD2 

0 . DOO , 0 . DOO , 0 . DOO , 

mll,ml2,ml3,m21, m22 ,m23 ,m31 , m32 ,m33 , ml ,m2 ,m3) 


ELSE 

* 


IF (I . EQ. 1)THEN 
* Mmc=Mms*Msc 


c 


* 

c 


* 


eq. (2.26) 

CALL COMB I (m 1 1 , m 1 2 , m 1 3 , m2 1 , m22 , m23 , m3 1 , m32 , m3 3 , ml , m2 , m3 , 

1 . DOO , 0 . DOO , 0 . DOO , 0 . DOO , CSPHIG , -SNPHIG , 0 . DOO , SNPHIG CSPHIG 

0. DOO, 0. DOO, 0. DOO, 

1. DOO, 0. DOO, 0. DOO, 0. DOO, CSQG,SNQG,0. DOO, -SNOG, CSQG, 

O.DOO, SG*SNQG, SG*CSQG) 

END IF 


Mpc=Mpm ,lc Mmc 


eqs. (2.25), (3.13) 

CALL COMBI (pll.pl2.pl3, p 21,p22,p23,p31, P 32,p33, pi, P 2,p3, 

CSD2 , 0 . DOO , -SND2 , 0 . DOO , 1 . DOO , 0 . DOO , SND2 , 0 . DOO , CSD2 
0. DOO, 0. DOO, 0. DOO, 


mll,ml2,ml3,m21,m22,m23,m31,m32,m33,ml ,m2,m3) 

END IF 


* DETERMINE MAIN CONTACT POINT 

* 


* 

* CALCULATE THETAG 

* 


c X ( 1) represents THETAG 

PAR (2) = (MG2-SNRT2) *CSPSIG 
IF(HG .EQ. 1 L ' ) THEN 
PAR (3)— SNQG*CSRT2*SNPSIG 
c step 1 in sec. 3.2 

X (1) =QG-BETA+RITAG 
ELSE 

PAR (3)=SNQG*CSRT2*SNPSIG 
c step 1 in sec. 3.2 

X (1) =- (QG-BETA+RITAG) 

END IF 

CALL NONLIN (PCN, 14,1,100,X,F,FI, 1 . D~5 , A2SP , IPVTP WORKP) 
THETAG=X (1) 

CSTHEG=DCOS (THETAG) 

SNTHEG=DS IN (THETAG) 

* 

* CALCULATE TAUG 


c eq. (2.38) 

IF(HG .EQ. ' L ' ) THEN 
TAUG=THETAG-QG+PHIG 
ELSE 

TAUG=THETAG+QG-PHIG 
END IF 

CSTAUG=DCOS (TAUG) 
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SNTAUG=DSIN (TAUG) 


* CALCULATE UG 

* 

c. eq. (2.43) 

IF (HG .EQ. ' L ' ) THEN 

UG=RG*CTPSIG*CSPSIG-SG* ( (MG2-SNRT2) *CSPSIG*SNTHEG-DSIN (QG-PHIG) * 
// CSRT2*SNPSIG) / (CSRT2*SNTAUG) 

ELSE 

UG=RG*CTPSIG*CSPSIG-SG* ( (MG2-SNRT2) *CSPSIG*SNTHEG+DSIN (QG-PHIG) * 
# CSRT2*SNPSIG) / (CSRT2*SNTAUG) 

END IF 

* 

* CALCULATE MAIN CONTACT POINT 

* 

c eq. ( 2 . 1) 

Bcx=RG*CTPSIG-UG*CSPSIG 

Bcy=UG*SNPSIG*SNTHEG 

Bcz=UG*SNPSIG*CSTHEG 

c e q . (2.2) 

Ncx=SNPSIG 

Ncy=CSPSIG*SNTHEG 

Ncz=CSPSIG*CSTHEG 

c eq . (2.9) 

EGIcx=0.D00 

EGIcy=CSTHEG 

EGIcz=-SNTHEG 

c eq . (3.13) 

CALL TRCOOR (Bpx , Bpy, Bpz, 

. pll,pl2,pl3,p21 ,p22,p23,p31,p32,p33,pl ,p2,p3, 

. Bcx,Bcy,Bcz) 

c eq. (3.16) 

CALL TRCOOR (Npx , Npy ,Npz, 

. pll,pl2,pl3,p21,p22,p23,p31,p32,p33,0. D00 , 0 . D00 , 0 . D00, 

. Ncx,Ncy,Ncz) 

c eq. (3.17) 

CALL TRCOOR (EGIpx , EGIpy, EGIpz, 

. pll,pl2,pl3,p21,p22,p23,p31,p32,p33, 0. D00, 0.D00, 0. D00, 

. EG lex , EGIcy , EGIcz) 

c fig. 18 & sec. 3.3 

Bf x=Bpx 
Bfy=Bpy 
Bfz=Bpz 
Nf x=Npx 
Nfy-Npy 
Nf z = Npz 
EGIf x=EGIpx 
EGIfy=EGIpy 
EGI f z=EGIpz 

* 

* CALCULATE PS IP 

* 

c eq. (3.83) 
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PSIP=DASIN(CSDl*Nfx-SNDl*Nfz) 

IF (I .EQ. 1)THEN 
PSIP=-PSIP+180.D00*CNST 
END IF 

CSPSIP=DCOS (PSIP) 

SNPSIP=DSIN (PSIP) 

* 

* CALCULATE TAUP 

* 

c eqs. (3 .84) - (3 .86) 

TAUP=DATAN2 (Nfy/CSPSIP, (Nf x-CSDl*SNPSIP) / (~SND1*CSPSIP) ) 
CSTAUP=DCOS (TAUP) 

SNTAUP=DSIN (TAUP) 

* 

* CALCULATE PRINCIPAL CURVATURES AND DIRECTIONS OF THE GEAR CUTTER 


c eq. (2. 10) 

KGI=-CTPSIG/UG 

c eq. (2.12) 

KGII=0.D00 

c the second principal direction is determined by rotating of 

c the first principal derection about unit normal by 90 degrees 


CALL ROTATE (EGI If x , EGI I f y , EGI I f z, EGIf x , EGIf y , EG I f z, RITAG , 
. Nfx.Nfy ,Nfz) 

* 

* CALCULATE W2G 

* 

c eqs. (3 . 18) - (3 . 20) 

IF (HG .EQ. ' L ' ) THEN 
W2fx=-SNPIT2 
WGfx— MG2*CSD2 
W2fy=0.D00 
WGfy=0.D00 
W2fz=CSPIT2 
WGfz*-MG2*SND2 
ELSE 

W2fx=SNPIT2 

WGfx*MG2*CSD2 

W2fy=0.D00 

WGf y=0 . D00 

W2fz=-CSPIT2 

WGfz-MG2*SND2 

END IF 

* 

W2Gf x=W2f x-WGfx 
W2Gf y=W2f y-WGf y 
W2Gfz=W2fz-WGfz 

Vc 

* CALCULATE VT2 , VTG, AND VT2G 

* 

c eq. (3.22) 

CALL CROSS (VT2f x , VT2fy , VT2fz, W2fx, W2fy, W2fz,Bf x , Bfy , Bf z) 

c eq. (3.21) 
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CALL CROSS (VTGf x , VTGfy , VTGf z, WGf x , WGfy , WGf z, Bfx, Bfy,Bfz) 

eq . (3.23) 

VT2Gf x=VT2f x-VTGf x 
VT2Gfy=VT2fy - VTGfy 
VT2Gfz=VT2fz-VTGfz 

CALCULATE V(2G)GI AND V(2G)GII 

eq. (3.24) 

CALL DOT (VGI , EGIf x , EGIfy , EGIf z, VT2Gf x, VT2Gfy ,.VT2Gf z) 

eq. (3.25) 

CALL DOT (VGI I , EGI I f x , EGI I f y , EGI I f z , VT2Gf x , VT2Gfy , VT2Gfz) 

CALCULATE Al3,A23,A33 
eq. (3.26) 

CALL DET (DETI , W2Gfx , W2Gfy,W2Gfz, Nfx , Nfy, Nfz, EGIf x , EGI fy, EGI fz) 
A13=-KGI*VGI-DETI 
eq. (3.27) 

CALL DET (DETII , W2Gf x , W2Gfy , W2Gfz, Nf x,Nfy ,Nf z.EGIIf x , EGI Ify , EGIIf z) 
A23=-KGII 5V VGII-DETII 
eq. (3.28) 

CALL DET (DET3 , Nf x , Nfy , Nf z, W2Gf x , W2Gf y , W2Gf z, VT2Gf x , VT2Gfy, VT2Gf z) 
CALL CROSS (Cx,Cy,Cz,W2fx,W2fy,W2fz, VTGf x, VTGfy, VTGfz) 

CALL CROSS (Dx , Dy , Dz , WGf x , WGfy , WGf z , VT2f x , VT2f y , VT2f z) 

CALL DOT (DET45 ,Nfx,Nfy,Nfz, Cx~Dx , Cy-Dy , Cz~Dz) 
A33=KGI*VGI*VGI+KGII ,v VGII*VGII-DET3-DETA5 

CALCULATE SIGMA 

eq. (3.29) 

P=A23*A23“Al3*Al3+ (KGI-KGII) *A33 
SIGDBL=DATAN (2 . D00*Al3*A23/P) 

SIGMA=0.5D00*SIGDBL 

CALCULATE K2I AND K2II 

eqs. (3 . 30) - (3 . 31) 

Tl=P/ (A33*DCOS (SIGDBL) ) 

T2=KGI+KGII- (A13*A13+A23*A23) /A33 
K2I= (T1+T2) /2 . D00 
K2II=(T2-T1)/2.D00 

CALCULATE E2I AND E2II 

description after eq. (3.29) 

CALL ROTATE (E2lf x, E2I f y, E2If z, EGIf x, EGIfy, EGIf z, -SIGMA, Nfx, Nfy, 

. Nfz) 

CALL ROTATE (E2I If x , E2I I f y , E2I If z, E2lf x , E2lfy , E2If z, RITAG, 

. Nfx, Nfy, Nfz) 

eq. (3. 44) 

TNETAG=DSIN (DLTA+SIGMA) /DCOS (DLTA+SIGMA) 



CALCULATE W2 


eq. (3.33) 

IF (HG .EQ. ' L ' ) THEN 
W2fx=-MM21*SNPIT2 
W2fy=0 . DOO 
W2fz=MM21*CSPIT2 
ELSE 

W2fx=MM21*SNPIT2 
W2fy=0.D00 
W2fz=-MM21*CSPIT2 
END IF 

CALCULATE W1 

eq. (3.32) 

IF (HG .EQ. ' L ' ) THEN 
Wlfx=-SNPIT1 
Wlfy=O.DOO 
Wlfz=-CSPIT1 
ELSE 

Wlfx=SNPITl 
Wlfy=O.DOO 
Wlf z=CSPITl 
END IF 

CALCULATE W12 

eq. (3.34) 

Wl2fx=Wlfx-W2fx 

W12fy=Wlfy-W2fy 

Wl2fz=Wlfz-W2fz 

CALCULATE VT2 

eq. (3.36) 

CALL CROSS (VT2fx , VT2f y , VT2f z, W2f x , W2f y , W2fz, Bf x , Bfy , Bfz) 
CALCULATE VT1 
eq. (3.35) 

CALL CROSS (VTlf x, VTlfy, VTlf z, Wlf x.Wlfy.Wlfz.Bfx, Bfy, Bfz) 

CALCULATE VT12 

eq. (3.37) 

VT12fx=VTlfx-VT2fx 

VTl2fy=VTlfy-VT2fy 

VTl2fz=VTlfz-VT2fz 


CALCULATE V2 


c 


CALL DOT (V2l,VTl2fx,VT12fy,VTl2fz,E2lfx, E2Ify, E2lfz) 

- eq. (3.39) 

CALL DOT(V2II,VT12fx,VT12fy,VT12fz,E2IIfx,E2IIfy,E2lIfz) 

* CALCULATE A31 

* 

c eq. (3.40) 

CALL DET (DET1 ,W12fx,Wl2fy,W12fz,Nfx,Nfy,Nfz,E2lfx,E2Ify,E2lfz) 
A31=-K2I ,V V2I-DET1 

c eq. (A. 33) 

Al3=A31 

* 

* CALCULATE A3 2 

jV 

c eq. (3.41) 

CALL DET(DET2,W12fx,Wl2fy,Wl2fz,Nfx,Nfy,Nfz,E2IIfx,E2IIfy,E2lIfz) 
A32=-K2 I I *V2 I I -DET2 

c eq. (A. 35) 

A23=A32 

* 

* CALCULATE A3 3 

* 

c eq. (3.42) 

CALL DET(DET3,Nfx,Nfy,Nfz,Wl2fx, Wl2fy, Wl2fz, VTl2fx,VT12fy, VTl2fz) 
CALL CROSS (Cx,Cy,Cz,Wlfx, Wlfy, Wlfz, VT2fx, VT2fy, VT2fz) 

CALL CROSS (Dx,Dy,Dz, W2fx, W2fy, W2fz, VTlfx, VTlfy, VTlfz) 

CALL DOT (DOT1 , Nf x ,Nf y , Nf z, Cx-Dx , Cy-Dy , Cz~Dz) 

CALL DET (DET4,Nf x ,Nfy , Nf z, W2fx ,W2fy,W2fz,Bfx,Bfy,Bfz) 

A33=K2 I*V2 I *V2 I +K2I I*V2 I I *V2 I I -DET3-DOT 1 +M2 1 PRM*DET4 

* 

* CALCULATE ETAP 

* 

c eq. (3.53) 

ETAP=DATAN ( ( (A33+A3 1 *V2 1) *TNETAG"A3 1 *V2I I) / (A33+A32* 

. (V2II-V2I*TNETAG))) 

TNETAP=DSIN (ETAP) /DCOS (ETAP) 

* 

* CALCULATE All, Al2, AND A22 

* 

N3= ( 1 . DOO+TNETAP*TNETAP) *A33 

c eq. (3.72) 

Nl= (A13*A13- (A23*TNETAP) **2) /N3 
c eq. (3.73) 

N2= (A23+A13 ,V TNETAP)*(A13+A23*TNETAP)/N3 

KS2=K2I+K2II 

G2=K2I-K2II 

c eqs . (3.74), (3.75) 

KS 1=KS2- ( (4 . D00*AXI A*AXIA-N 1 *N 1 -N2 ,V N2) * ( 1 . DOO+TNETAP*TNETAP) / 

. (-2 . D00*AXIA* (1 . DOO+TNETAP*TNETAP) +N1* (TNETAP*TNETAP-1 . D00) 

. -2.D00*N2*TNETAP)) 


c eqs. (3.66), (3.69) & description after eq. (3.60) 

Al 1=TNETAP*TNETAP/ (1 . DOO+TNETAP*TNETAP) * (KS2-KS1) +N1 
c eqs. (3.67), (3.70) & description after eq. (3.60) 
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A12=-TNETAP/ (1 . D00+ TNETAP*TNETAP) * (KS2-KS1) +N2 

c eqs. (3.68), (3.71) & description after eq. (3.60) 

A22= 1 . D00/ ( 1 . DOO+TNETAP*TNETAP) * (KS2-KS1) ~N1 

c eq. (A. 32) 

A21=A12 


* CALCULATE SIGMA (12) 

* 

c eq. (3.77) 

SIGDBL=DATAN (2 . D00*A12/ (K2I-K2I I-Al 1+A22) ) 

SIGM12=.5D00*SIGDBL 

* 

* CALCULATE KlI AND Kill 

* 

c eq. (3.78) 

Gl=2.D00*Al2/DSIN (SIGDBL) 

c eq. (3.79) 

KlI=.5D00 ,,f (KS1+G1) 

KlI 1= . 5D00* (KS1-G1) 

* CALCULATE Ell AND El I I 

Vc 

c similar to description after eq. (3.29) 

CALL R0TATE(ElIfx,ElIfy,ElIfz,E2lfx,E2lfy,E2lfz,-SIGM12,Nfx,Nfy, 
. Nf z) 

CALL ROTATE (ElIIfx,ElIIfy,ElIIfz, El If x, El If y , El If z, RITAG, 

. Nfx,Nfy,Nfz) 

it 

* PINION 


* CALCULATE PRINCIPAL DIRECTIONS OF THE PINION CUTTER 

y. 

c eq. (3.92) 

IF (HG .EQ. 'L')THEN 
EPIf x"SNDl*SNTAUP 
EPIfy=CSTAUP 
EPIfz=CSDl’ v SNTAUP 
ELSE 

EP I f x=-SNDl *SNTAUP 
EPIfy=-CSTAUP 
EPIf z=~CSDl *SNTAUP 
END IF 

IF(DACOS(EGIfx*EPIfx+EGIfy*EPIfy+EGIfz*EPIfz)/CNST. .GT. A5.D00) 
.THEN 

EPIfx=-EPIfx 
EPIfy=-EPIfy 
EPIfz=-EPIfz 
END IF 

CALL ROTATE (EPI If x , EPIIf y , EPIIf z, EPIf x, EPIf y, EPIf z, RITAG, 

. Nfx,Nfy,Nfz) 

?v 

* CALCULATE THE ANGLE BETWEEN PRINCIPAL DIRECTIONS OF PINION AND CUTTER 
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it 

cross product of eli and epi 

SNSIGM-(ElIfy*EPIfz-ElIfz*EPIfy)/Nfx 

. dot product of eli and epi 

CSSIGM=Ellfx*EPIfx+ElIfy*EPIfy + ElIf z*EPIf z 

CS2SIG=2 . DOO*CSSIGM*CSSIGM~ 1 . D00 

TN2SIG=2.D00*SNSIGM*CSSIGM/CS2SIG 

* CALCULATE PRINCIPAL CURVATURES OF PINION CUTTER 

c eq. (2.12) 

KPII=0.D00 

C w!iimin/«H*SKSI<M*SMSI®«in*CSSIOH*CSSIGK) 

it 

* CALCULATE All, Al2, AND A22 

it 

C j^ll.KPI-KlI*CSSIGH*CSSIGM-KlII*SNSIGM*SNSIGM 

c eq. (A. 32) 

Al2= (K1I-K1II)*SNSIGM*CSSIGM 

C A 22=KPII-KlI , ''SNSIGM’ t SNSIGM-KlII’ , 'CSSIGM ,t CSSIGH 

* 

* CALCULATE UP 

it 

c eq. (3.95) 

UP=1.D00/(KPI*SNPSIP/CSPSIP) 

* 

* CALCULATE RP 

it 

C eq. (3.99) 

Bmx=-Bfx*CSDl+Bfz*SNDl 

eq . (3.100) 

RP= (Bmx+UP’' ! CSPSIP) *SNPSIP/CSPSIP 

it 

* CALCULATE MPl 

C C 22 --w)*CT > nfr-N£^EPII£y)*SNPITU(»£y*EPIIE»-NEx*EPnE y )*CSPITl 

IF (HG .EQ. 1 R ' ) THEN 
C11=-C11 
C12=-C12 
C22=-C22 
END IF 

/o 119) 

T4 =< (Bfy*CSRTl) / (EPIIf x*CSDl-EPIIfz*SNDl) 

IF (HG .EQ. 1 R 1 ) THEN 
T4=-T4 
END IF 
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C — — eq. (3.120) 

T1=-C11/KPI 

c I 2= (A1 J* K ^J* T4+A11 * C22 ' A12 * C12) / (A12*KPI) 

• vJ • 122) 

Ull=Tl*EPIfx 

U12=T2*EPIfx+T4*EPIIfx 

U21=Tl*EPIfy 
U22=T2*EPIfy+T4*EPIIfy 
U31=Tl*EPIf z 
U32=T2*EPIfz+T4*EPIIfz 
c eq. (3.124) 

c ! 1 ^ 21 )3f^ ) SD1 * U21 ' ,Nf **SNDl-NfynuU*SNDl*U31»CSDl) 

VI— VI 
V2=-V2 
V3=-V3 
END IF 

c eq. (3.132) 

Hll=-U21 ,t CSPITl+SNDl*(Bfz*SNPITl-Bfx*CSPTT1'> 

c eq. (3.134) ' 

H21=Ull*CSPITl-U31*SNPITl+Bfy*SNRTl 

c eq . (3 . 136) 

c H31=U21*SNPIT1+CSD1* (Bfz*SNPITl-Bfx*CSPITl) 

H12= (Bfz*SNPITl-Bfx*CSPITl-U22) *CSPIT1 
c eq. (3.135) 

H22— (Bfy-U12*CSPIT1+U32*SNPIT1) 
c eq. (3.137) 

H32=- (Bfz :t SNPITl-Bfx*CSPITl-U22) *SNPIT 1 
IF(HG .EQ. 1 R 1 ) THEN 

H11-U21 X CSPIT1+SND1* (Bf z*SNPITl~Bf x*CSPITl) 

H21 LJ1 l 5 ''CSPITl+U31 5V SNPITl+Bfy ,v SNRTl 

H31=-U 21,f SN piT l+CSDl*( B fz*SN p ITl-Bfx*CSPITl) 

H12- (Bfz ,t SNPITl-Bfx*CSPITl+U22) *CSPIT1 
H22— (Bfy+U12*CSPIT1-U32*SNPIT1) 
H32=-(Bfz*SNPITl-Bfx*CSPITl+U22)*SNPITl 

END IF 1 

c eq. (3.139) 

Fl=Nfx*Hll+Nfy*H21+Nfz*H31 
c eq. (3.140) 

F2=Nfx*H12+Nfy*H22+Nfz*H32 

c eq. (3.145) 

Y2=A12*(2.D00*KPI*T1*T2-V2-F1) 

Sl--W« I * I2 * T2 * KPII ' , "‘ T4 - V3 - r2) -( KPI * T2t C12)*<KPII*I^C22) 

,c CALCULATE EM AND LM 
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C eq. (3.122) 

VTIPf x=Ul 1 ,V MP1+U12 
VTlPfy=U21*MPl+U22 
VTlPfz=U31 ,v MPl+U32 

c eq. (3.111) 

IF (HG .EQ. ' L ' ) THEN 
EM= (Bfy*CSPITl-VTlPfx) / (MP r' c SNDl) +Bf y 

LM=(Bfx*CSPITl-Bfz*SNPITl+VTlPfy)/MPl+Bfx*SNDl+Bfz*CSDl 

ELSE 

EM= (-Bfy*CSPITl-VTlPfx)/(MPl y, SNDl)-Bfy 

LM=(Bfx*CSPITl-Bfz*SNPITl-VTlPfy) /MPl+Bf x*SNDl+Bf z*CSDl 
END IF 

* 

* CALCULATE SP AND QP 

iV 

C eqs. (3.150) , (3.151) 

IF (HG .EQ. ' L ' ) THEN 
Zl— Bfy+EM-UP*SNPSIP*SNTAUP 
ELSE 

Zl=Bfy+EM+UP*SNPSIP*SNTAUP 
END IF 

Z2=Bfx*SNDl+Bfz*CSDl-LM-UP*SNPSIP*CSTAUP 

SP=DSQRT(Z1*Z1+Z2*Z2) 

QP=DATAN(Z1/Z2) 

IF (HG .EQ. ' L ' ) THEN 
THETAP=TAUP~QP 
ELSE 

THETAP=TAUP+QP 
END IF 

it 

* CONVERT RADIAN TO DEGREE 

* 

PSIGDG=PSIG/CNST 

PSIPDG=PSIP/CNST 

TAUGDG=TAUG/CNST 

TAUPDG=TAUP/CNST 

QGDG=QG/CNST 

QPDG=QP/CNST 

THEGDG= THETAG/CNST 

THEPDG= THET AP / CNS T 

PHIGDG=PHIGO/CNST 

* 

* OUTPUT 

* 

WRITE (72 , 1 0000) PS IGDG , PS I PDG , RG , RP , TAUGDG , TAUPDG , SG , SP , QGDG , QPDG , 
. MG2 , MP 1 , EM , LM , UG , UP , THEGDG , THEPDG , PH I GDG 
10000 FORMAT (IX, 'PSIGDG = ' , G20 . 12 , 12X , ' PSIPDG -'.G20.12,/ 

, IX, 1 RG = ' , G20 . 12 , 12X, ' RP =',G20.12,/ 

, IX, 'TAUGDG =' ,G20. 12, 12X, 'TAUPDG =',G20.12,/ 

, IX, ' SG =' ,G20. 12, 12X, 1 SP =',G20.12,/ 

, IX, 'QGDG =’ ,G20.12,12X, ' QPDG =',G20.12,/ 

, IX, ' MG2 =' ,G20.12,12X, 'MPl =',G20.12,/ 

, IX, 'EM = ' ,G20. 12, 12X, ' LM =',G20.12,/ 
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, IX, ' UG =' ,G20. 12, 12X, 'UP =',G20.12,/ 

. , IX, ' THETAGDG = ' ,G20. 12, 12X, 'THETAPDG =',G20.12,/ 

, IX, ' PHIGODG =' ,G20.12,12X,/) 

?v 

* TCA 

* 


IF (I .EQ. 1)THEN 

TPAR(1)=RG*CSPSIG/SNPSIG*CSPSIG 

TPAR (2) = (MG2-SNRT2) *CSPSIG 

TPAR (3) =CSRT2*SNPSIG 

TPAR (4) =RG*CSPSIG/SNPSIG 

TPAR (5) =CSD2*SNPSIG 

TPAR (6) =SND2*CSPSIG 

TPAR (7) =SND2*SNPSIG 

TPAR (8) =CSD2*CSPSIG 

TPAR (9) =RP*CSPSIP/SNPSIP*CSPSIP 

TPAR (10) = (MPl-SNRTl) *CSPSIP 

TPAR (11) =CSRTl*SNPSIP 

TPAR (12) =SNRT1 ,C CSPSIP 

TPAR (13) =RP*CSPSIP/SNPSIP 

TPAR(14)-CSD1*SNPSIP 

TPAR (15) =SND1*CSPSIP 

TPAR (16) =SND1*SNPSIP 

TPAR (17) =CSD1*CSPSIP 

TPAR (18) =LM*SND1 

TPAR (19) =LM*CSD1 

PHIP=0.D00 
PHI21=0.D00 
PHI 1 1=0. D00 
CSPH1 1=DC0S (PHI 1 1) 

SNPH1 1=DSIN (PHI 1 1) 

TX (1) =PHIP 
TX (2) =THETAP 
TX (3) =PHI21 
TX (4) =PHIGO 
TX (5) =THETAG 

CALL NONLIN (TCN , 14 , 5 , 100 , TX, TF , TF1 , 1 . D~5 , AZS , IPVT , WORK) 
PHIP0=TX (1) 

THEP0=TX (2) 

PHI210=TX (3) 

PHIG0=TX (4) 

THEG0=TX (5) 

TX (1) =PHIP0 
TX(2) =THEP0 
TX (3) =PHI210 
TX(4) =PHIG0 
TX(5) =THEG0 

D1HI11=18.D00/36.D00*CNST 
DO 100 I J=1 ,60 
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CSPHll=DCOS (PHIll) 

SNPHl 1=DSIN (PHI 1 1) 

CALL NONLIN (TCN , 1 4 , 5 , 1 00 , TX , TF , TF 1 , 1 . D~5 , AZS , IPVT , WORK) 

PHIP=TX (1) 

THETAP=TX (2) 

PHI21=TX (3) 

PHIG=TX (4) 

THETAG=TX (5) 

ERROR= ( (PH I 2 1*36 . D02-PHI210*36 . D02) -PHI 1 l*36.D02*TNl/TN2) /CNST 

* 

CALL PRING2(KS2,G2,E2lfx,E2Ify,E2lfz,E2lIfx,E2IIfy,E2IIfz) 

CALL PRINP1 (KS1 ,Gl,ElIfx,ElIfy,ElIfz,ElIIfx,ElIIfy,ElIIfz) 

CALL SIGAN2 (E2If x , E2lfy,E2Ifz,E2lIfx,E2IIfy,E2llfz,ElIfx,ElIfy, 
ElIfz,CS2SIG,SN2SIG,SIGM12) 

CALL EULER (KS2,G2,KS1 , G1 , CS2SIG , SN2SIG, IEU) 

IF (IEU .EQ. 1)THEN 
WRITE (72,*) ’THERE IS INTERFERENCE' 

GO TO 88888 
END IF 

CALL ELL I PS (KS2 , G2 , KS1 , G1 , CS2SIG, SN2SIG, DEF , ALFAl , 

AXISL , AXISS ,ElIfx,ElIfy,ElIfz) 

* 

CALL PF (B2px , B2py ,B2pz,B2fx,B2fy,B2fz) 

* 

* XBf, YBf, and ZBf is the direction of the long axis of the ellipse 

Vc 

CALL PF (XBp , YBp , ZBp , XBf , YBf , ZBf ) 

ELBlpx=B2px+XBp 

ELBlpz=B2pz + ZBp 

ELB2px=B2px“XBp 

ELB2pz=B2pz~ZBp 

v? 

IF (I .EQ. 1)THEN 

WRITE (79 , 9000) I J , PHI 1 1 /CNST , I J , ERROR 
WRITE (78, 8000) IJ,B2pz, IJ,B2px 
WRITE (77, 7000) ELBlpz,ELBlpx,ELB2pz,ELB2px 
ELSE 

WRITE (89, 9000) IJ, PH 111/ CNST, I J, ERROR 
WRITE (88, 8000) I J , B2pz, I J , B2px 
WRITE (87 , 7000) ELBlpz, ELBlpx , ELB2pz, ELB2px 
END IF 

* 

PHI 1 1=PHI 1 1+D1HI 1 1 

* 

100 CONTINUE 

it 

if 

if 

PHI11=O.DOO 
CSPH1 l=DCOS (PHIll) 

SNPH11=DSIN (PHI 11) 
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TX(1)=PHIP0 
TX (2) =THEPO 
TX (3) =PHI210 
TX (4) =PHIGO 
TX (5) =THEGO 

D1HI11=18.D00/36.D00*CNST 

X 

DO 200 IJ-1,60 
CSPHl l=DCOS (PHI 11) 

SNPH1 1=DSIN (PHI 1 1) 

CALL NONLIN (TCN , 1 4 , 5 , 1 00 , TX , TF , TF 1 , 1 . D~5 , AZS , IPVT , WORK) 

PHIP=TX (1) 

THETAP=TX (2) 

PHI21=TX (3) 

PHIG=TX (4) 

THETAG=TX (5) 

ERROR= ( (PHI2 1*36 . D02-PHI210*36 . D02) -PHI 11*36 . D02*TN1/TN2) /CNST 

X 

CALL PRING2 (KS2 , G2,E2lfx,E2lfy,E2lfz,E2lIfx,E2IIfy,E2lIfz) 

CALL PRINP1 (KS1 ,G1 , El If x , El Ify , El If z, El I If x , El I If y , El 1 1 f z) 

CALL SIGAN2 (E2Ifx,E2lfy,E2lfz,E2lIfx,E2IIfy,E2IIfz,ElIfx,ElIfy, 
ElIfz,CS2SIG,SN2SIG,SIGMl2) 

X 

CALL EULER (KS2 , G2 , KS1 , G1 , CS2SIG, SN2SIG, IEU) 

IF (IEU .EQ. 1)THEN 
WRITE (72,*) 'THERE IS INTERFERENCE' 

GO TO 88888 
END IF 

* 

CALL ELL I PS (KS2 , G2, KS1 , G1 , CS2SIG, SN2SIG.DEF, ALFAl , 

AXISL , AXISS , El If x , El If y , Ellfz) 

* 

CALL PF (B2px , B2py , B2pz, B2f x , B2fy , B2f z) 

* 

* XBf, YBf, and ZBf is the direction of the long axis of the ellipse 

Vc 

CALL PF (XBp , YBp , ZBp , XBf , YBf ,ZBf ) 

ELBlpx=B2px+XBp 

ELBlpz=B2pz + ZBp 

ELB2px=B2px~XBp 

ELB2pz“B2pz-ZBp 

* 

IF (I .EQ. 1)THEN 

WRITE (79 , 9001) I J, PHI 1 1/CNST , I J , ERROR 
WRITE (78, 8001) I J,B2pz, I J,B2px 
WRITE (77, 7000) ELBlpz,ELBlpx,ELB2pz,ELB2px 
ELSE 

WRITE (89 , 9001) I J , PHI 1 1/CNST , I J , ERROR 
WRITE (88,8001) I J , B2pz, I J , B2px 
WRITE (87, 7000) ELBlpz, ELBlpx , ELB2pz, ELB2px 
END IF 

sV 

PHI 1 1=PHI 1 1-DlHI 1 1 


146 



200 CONTINUE 

* 

END IF 

* 

99999 CONTINUE 
88888 CONTINUE 

7000 FORMAT (6X, ' EX (1) = ' , F9 . 6 , / , 6X , 1 EY (1) = ' , F9 . 6 , / , 

6X, ’ EX (2) = ' , F9 . 6 , / , 6X, ' EY (2) = ’ , F9 . 6 , / , 

6X , 1 CALL CURVE (EX , EY , 2 , 0) 1 ) 

8000 FORMAT (6X, 'XO ( ' , 12, ’)= ' , F9 . 6 , / ,6X, ' YO ( ' , 12 , ' )= ' , F9.6) 

8001 FORMAT (6X, ' Xl ( ' , 12, ’ ) = ' , F9 . 6 , / ,6X, ' Y1 ( ' , 12, ' )= ' , F9.6) 

9000 FORMAT (6X, ' XO ( ' , 12 , ' ) = ' , F7. 3, / ,6X, ' YO ( ' , 12, ' ) , F8. 3) 

9001 FORMAT (6X, ' Xl ( ' , 12 , ' ) = ' , F7 . 3 , / , 6X, ' Yl ( 1 , 12, ' ) = ' , F8 . 3) 

END 

* for THE DETERMINATION OF MEAN CONTACT POINT 

it 

SUBROUTINE PCN(X,F,NE) 

IMPLICIT REAL*8(A-H,K,M-Z) 

CHARACTER*8 HG 
INTEGER NE 

REAL*8 X (NE) , F (NE) , PAR (6) 

COMMON/P 1 /PAR 
COMMON/AO/HG 

COMMON/ Al/p Il,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3 
COMMON/A2/SG, CSRT2 , QG, SNPSIG, CSPSIG 
COMMON/ A3 /TND1 , TND2 , R I TAG 
THETAG=X (1) 

CSTHEG=DCOS (THETAG) 

SNTHEG=DSIN (THETAG) 

IF (HG .EQ. ' L ' ) THEN 

UG=PAR (1) -SG* (PAR (2) *SNTHEG+PAR (3) ) / (CSRT2*DSIN (THETAG-QG) ) 
ELSE 

UG=PAR (1) -SG* (PAR (2) *SNTHEG+PAR (3) ) / (CSRT2*DSIN (THETAG+QG) ) 
END IF 

Bcx=PAR (4) -UG*CSPSIG 
Bcy=UG*SNPSIG*SNTHEG 
Bcz=UG*SNPSIG*CSTHEG 
CALL TRCOOR (Bpx , Bpy , Bpz, 

. pll,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3, 

. Bcx,Bcy,Bcz) 

XM=Bpz* (TND1-TND2) /2 . D00 
F(l)=Bpx-XM 

END 

y? 

* FOR THE DETERMINATION OF COORDINATES AND NORMALS OF CONTACT POINTS 

* 

SUBROUTINE TCN (TX , TF , NE) 

IMPLICIT REAL*8(A-H,K,M~Z) 

INTEGER NE 
CHARACTER *8 HG 
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REAL *8 TX (NE) , TF (NE) , TPAR (19) , LM 

COMMON/AO/HG 

COMMON/ T 1 / TPAR 

COMMON/A2/SG,CSRT2,QG,SNPSIG,CSPSIG 
COMMON /A4/CSD2 , SND2 , CSPIT2 , SNPIT2 

COMMON/B1/CSPH1 1 , SNPH1 1 , SP , EM, LM, CSRT1 , CSDl , SND1 , CSPSIP , SNPSIP 
COMMON/B2/CSPIT1 , SNPITl ,MP1 ,MG2,QP 

COMMON/B3/B2fx,B2fy,B2fz 
COMMON/B4/CSPH2, SNPH2,CSPH21 , SNPH21 
COMMON/C1/UG, CSTAUG, SNTAUG 
COMMON/C2/N2fx ,N2fy ,N2fz 
COMMON/DI /UP , CSTAUP , SNTAUP 
COMMON/ Fl/PHIGO 
common/gi/da.dv 

PHIP=TX (1) 

THETAP=TX (2) 

PHI21=TX(3) 

PHIG=TX (4) 

THETAG=TX (5) 

CSPHIP=DCOS (PHIP) 

SNPHIP=DSIN (PHIP) 

CSTHEP=DCOS (THETAP) 

SNTHEP=DSIN (THETAP) 

CSPH21=DCOS (PHI21) 

SNPH21=DSIN (PHI21) 

CSPHIG=DCOS (PHIG) 

SNPHIG-DSIN(PHIG) 

CSTHEG=DCOS (THETAG) 

SNTHEG=DSIN (THETAG) 

PHI2= (PHIG-PHIGO) /MG2 

PHI1=PHIP/MP1 

CSPH2=DCOS(PHI2) 

SNPH2=DSIN (PHI2) 

CSPHl=DCOS (PHI1) 

SNPH1=DSIN (PHI 1) 

IF (HG .EQ, ' L ' ) THEN 
TAUP=THETAP+QP-PHI.P 
ELSE 

TAUP=THETAP~QP+PHIP 
END IF 

CSTAUP=DCOS (TAUP) 

SNTAUP=DSIN (TAUP) 

IF (HG .EQ. ' L 1 ) THEN 
TAUG=THETAG-QG+PHIG 
ELSE 

TAUG=THETAG+QG-PHIG 
END IF 

CSTAUG=DCOS (TAUG) 

SNTAUG=DSIN (TAUG) 

CSQPHP=DCOS (QP-PHIP) 

SNQPHP=DSIN (QP-PHIP) 

CSQPHG=DCOS (QG-PHIG) 

SNQPHG=DSIN (QG-PHIG) 
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* GEAR 

* 

* SURFACE EQUATIONS 

* 

IF (HG .EQ. ' L ' ) THEN 

UG=TPAR (1) -SG* (TPAR (2) *SNTHEG-SNQPHG*TPAR (3) ) / (CSRT2*SNTAUG) 
B2py=UG*SNPSIG*SNTAUG-SG*SNQPHG 
ELSE 

UG=TPAR (1) -SG* (TPAR (2) *SNTHEG+SNQPHG*TPAR (3) ) / (CSRT2*SNTAUG) 
B2py=UG*SNPSIG*SNTAUG+SG*SNQPHG 
END IF 

B2px-CSD2* (TPAR (A) -UG*CSPSIG) -SND2* (UG*SNPSIG*CSTAUG+SG*CSQPHG) 
B2pz=SND2* (TPAR (A) -UG*CSPSIG) +CSD2* (UG*SNPSIG*CSTAUG+SG*CSQPHG) 
N2px=TPAR (5) -TPAR (6) *CSTAUG 
N2py=CSPSIG*SNTAUG 
N2pz=TPAR (7) +TPAR (8) *CSTAUG 

* 

* [Mwp] = [Mwa] [Map] 

it 

IF (HG .EQ. ' L ' ) THEN 

CALL COMB I (wp 1 1 , wp 1 2 , wp 1 3 , wp2 1 , wp22 , wp23 , wp3 1 , wp32 , wp33 , 

. wpl,wp2,wp3, 

. CSPH2 , SNPH2 , 0 . D00 , -SNPH2 , CSPH2 , 0 . D00 , 0 . DOO , 0 . DOO , 1 . DOO , 

. 0. DOO, 0. DOO, 0. DOO, 

. CSP IT2 , 0 . DOO , SNP IT2 , 0 . DOO , 1 . DOO , 0 . DOO , -SNP IT2 , 0 . DOO , CSP I T2 , 

. 0 . DOO , 0 . DOO , 0 . DOO) 

ELSE 

CALL COMB I (wp 1 1 , wp 1 2 , wp 1 3 , wp2 1 , wp22 , wp23 , wp3 1 , wp32 , wp33 , 

. wpl,wp2,wp3, 

. CSPH2 , -SNPH2 , 0 . DOO , SNPH2 , CSPH2 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

. 0 . DOO , 0 . DOO , 0 . DOO , 

. CSPIT2 , 0 . DOO , SNPIT2 , 0 . DOO, 1 . DOO, 0 . D00,~SNPIT2, O.DOO,CSPIT2 , 

. 0. DOO, 0. DOO, 0. DOO) 

END IF 

CALL TRCOOR (B2wx , B2wy , B2wz, 

. wpl 1 , wpl2 , wpl3 , wp21 , wp22 , wp23 , wp3 1 , wp32 , wp33 , wpl , wp2 , wp3 , 

. B2px , B2py , B2pz) 

CALL TRCOOR (N2wx,N2wy,N2wz, 

. wpl 1 , wp 12 , wp 1 3 , wp2 1 , wp22 , wp23 , wp3 1 , wp32 , wp33 , 0 . DOO , 0 . DOO , 0 . DOO , 

. N2px , N2py , N2pz) 

* 

* * * * * * * * * Vc * * * Vc * * * * * * * * * * * * * * * * * * * * * * »V * * * * * * * * * * * * * * * * * * * * * * rt * * * * * * * * * * * 
Vc 

* [Mfw] = [Mfa] [Maw] 

* 

fall=CSPIT2 
fal2=0. DOO 
fal3=-SNPIT2 
fa21=0.D00 
£a22=l . DOO 
f a23=0. DOO 
fa31=SNPIT2 
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fa32=0.D00 
fa33=CSPIT2 
fal=0. dOO 
f a2=0. dOO 
f a3=0. dOO 

IF (HG .EQ. ' L ' ) THEN 

CALL COMBI (fwl 1 , fwl2 , fwl3 , f w2 1 , fw22 , fw23 , fw31 , f w32 , f w33 , 

. fwl,fw2,fw3, 

. CSPIT2,O.DOO,-SNPIT2,O.DOO, 1 .D00,0.D00,SNPIT2,0.D00,CSPIT2, 

. 0 . DOO , 0 . DOO , 0 . DOO , 

. CSPH2 1 , -SNPH2 1 , 0 . DOO , SNPH2 1 , CSPH21 , 0. DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

. 0. DOO, 0. DOO, 0. DOO) 

ELSE 

CALL COMBI (fwll,fwl2,fwl3, f w21 , f w22 , fw23 , f w31 , fw32 , f w33 , 

. fwl,fw2,fw3, 

. CSPIT2,0.D00,-SNPIT2,0.D00, l.D00,0.D00,SNPIT2,0.D00,CSPIT2, 

. 0, DOO, 0. DOO, 0. DOO, 

. CSPH21 , SNPH2 1 , 0 . DOO , -SNPH2 1 , CSPH2 1 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

. 0. DOO, 0. DOO, 0. DOO) 

END IF 

CALL TRCOOR (B2fx,B2fy,B2fz, 

. fwl 1 , f wl2 , f wl3 , fw21 , f w22 , f w23 , f w3 1 , f w32 , f w33 , fwl , f w2 , f w3 , 

. B2wx , B2wy , B2wz) 

CALL TRCOOR (N2fx,N2fy,N2fz, 

. fwl 1 , fwl 2 , fwl 3 , fw21 , fw22 , fw23 , fw31 , fw32, fw33, O.DOO, O.DOO , O.DOO, 
. N2wx,N2wy,N2wz) 

* 

* PINION 

* 

* SURFACE EQUATIONS 

* 

IF (HG .EQ. ' L ' ) THEN 

UP=TPAR (9) - (SP* (TPAR ( 1 0) *SNTHEP+SNQPHP*TPAR (1 1) ) -EM* (TPAR (11)+ 
TPAR (12) *CSTAUP) -LM*TPAR (12) *SNTAUP) / (CSRT1*SNTAUP) 
Blpy=UP*SNPSIP*SNTAUP+SP*SNQPHP-EM 
ELSE 

UP=TPAR (9) - (SP* (TPAR (10) *SNTHEP~SNQPHP*TPAR (11)) +EM* (TPAR (11)+ 
TPAR (12) *CSTAUP) -LM*TPAR (12) *SNTAUP) / (CSRT1*SNTAUP) 
Blpy=UP*SNPSIP*SNTAUP-SP*SNQPHP+EM 
END IF 

Blpx-CSDl*(TPAR(13)-UP*CSPSIP)-SNDl*(UP*SNPSIP*CSTAUP+SP* 

CSQPHP) ~LM*SND1 

Blpz=SNDl* (TPAR (13) ~UP*CSPSIP) +CSD1* (UP*SNPSIP*CSTAUP+SP* 

CSQPHP) +LM*CSD1 

Nlpx— (TPAR (1 A) -TPAR (15) *CSTAUP) 

Nlpy=-CSPSIP*SNTAUP 

Nlpz =_ (TPAR (16) +TPAR (17) *CSTAUP) 

* 

* [Mwp] = [Mwa] [Map] 

* 

IF (HG .EQ. 1 L 1 ) THEN 

CALL COMBI (wpl 1 , wpl2 , wpl3 , wp21 , wp22 , wp23 , wp31 , wp32 , wp33 , 

. wpl,wp2,wp3, 
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CSPH 1 , -SNPH1 , 0 . DOO , SNPH 1 , CSPHl , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

CSP I T i!o D DOO° SNP I T 1 , 0 . DOO , 1 . DOO , 0 . DOO , - SNP I T 1 , 0 . DOO , CSP I T 1 , 

0. DOO, 0. DOO, 0. DOO) 

CALL COMBI (wpll,wpl2,wpl3,wp21,wp22,wp23,wp31,wp32,wp33, 

CSPHl, SNPHl!o’.DOo!-SNPHl, CSPHl, 0. DOO, 0. DOO, 0. DOO, 1. DOO, 

0 CSm;!S?s“?h,O.DOO,l.DOO,O.DOO,-SNPITl,O.DOO,CSPITl, 

0. DOO, 0. DOO, 0. DOO) 

END IF 

TALL TRCOOR (Blwx , Blwy , Blwz, , 

. Wpll,wpl2,wpl3,wp21,wp22,wp23,wp31,wp32,wp33,wpl,wp2,wp3, 

. Blpx,Blpy,Blpz) 

, Nlpx ,Nlpy,Nlpz) 

;„,...*.».....**.*♦**♦*******»****•♦*”********************** 

* [Mpw] = [Mpa] [Maw] 

* 

X CALL COMB I (pwl l”pwl2 , pw 1 3 , pw2 1 , pw22 , pw23 , pw3 1 , pw32 , pw33 , 

CSPITl , (LDOO , -SNPITl , 0 .DOO, 1 .DOO, 0. DOO , SNP IT 1 ,O.DOO,CSPITl , 

; CSPHl i°SNPH 11 (hDOO , -SNPH1 1 , CSPHl 1 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

’ o. DOO, 0. DOO, 0. DOO) 

CALL COMBI (pwl 1 , pwl 2 , pwl 3 , pw21 ,pw22 , pw23, pw31 , pw32,pw33 , 

: Csim!o P DOO, -SNPITl, 0. DOO, 1. DOO, 0 . DOO, SNPIT! ,0.000, CSPITl , 

; CSPHl - SNPH 1 iToOO , SNPH 1 1 , CSPH 1 1 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

! 0. DOO, 6. DOO, 0. DOO) 

END IF 

CALL TRCOOR (Blpx,Blpy,Blpz, . - nu <i 

pwl 1 , P wl2 , pwl3 , pw2 1 , pw22 , pw23 , pw3 1 , pw32 , pw33 , pwl , pw2, pw3 , 

. Blwx, Blwy, Blwz) 

. ^pwl l[pwl2tpwl3*pw21^ pw22?pw23 , pw3 1 , pw32 , pw33 , 0 . DOO , 0 . DOO , 0 , DOO , 

. Nlwx,Nlwy,Nlwz) 

Blfx=-Blpx+DA*SNPIT1 

Blfy=-Blpy + DV 

Blfz=Blpz+DA*CSPITl 

Nlfx=-Nlpx 

Nlfy=-Nlpy 

Nl£z=Nlpz 

TF (1) =B2f x-Blfx 

TF(2)=B2fy-Blfy 

TF (3) =B2f z - Blfz 
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TF (4) =N2f x-Nlfx 
TF(5)=N2fy-Nlfy 
END 

FOR THE DETERMINATION OF GEAR PRINCIPAL CURVATURES AND DIRECTIONS 


SUBROUTINE PRING2 (KS2 , G2 , E2lf x 
IMPLICIT REAL*8(A-H,K,M-Z) 
CHARACTER *8 HG 
COMMON/AO/HG 


E2lfy,E2lfz,E2lIfx ) E2lIfy,E2lIfz) 


COMMON/A2/SG , CSRT2 , QG , SNPS IG , CSPS IG 
COMMON /A3 /TND 1 , TND2 , RITAG 

COMMON/ A4/ CSD2 , SND2 , CSP I T2 , SNP I T2 
COMMON/B2/CSPIT1 , SNPIT1 ,MP1 ,MG2 OP 
COMMON/B3/Bfx, Bfy, Bfz 
COMMON/BA/CSPH2 , SNPH2 , CSPH21 , SNPH21 
COMMON/ C 1 /UG , CSTAUG , SNTAUG 
COMMON/ C2/Nfx,Nfy,Nfz 
KGI=-CSPSIG/ (UG*SNPSIG) 

KGII=0. DOO 


EGIf x=SND2 ,,c SNTAUG 

EGIfy=CSTAUG 

EGIfz=-CSD2*SNTAUG 

CALL ROTATE (EGIIf x.EGIIfy , EGIIfz, EGIf* 
Nfx.Nfy.Nfz) 


EGIfy.EGIfz, RITAG, 


* CALCULATE W2G 

>v 

IF (HG .EQ. 1 L ' ) THEN 
W2fx=-SNPIT2 
WGfx=-MG2*CSD2 
W2f y=0 . DOO 
WGfy=0. DOO 
W2fz=CSPIT2 
WGfz— MG2*SND2 
ELSE 

W2f x=SNPIT2 
WGf x=MG2 ,t CSD2 
W2fy=0.D00 
WGf y=0 . DOO 
W2fz=-CSPIT2 
WGf z=MG2*SND2 
END IF 

W2Gf x=W2f x-WGfx 
W2Gfy=W2fy-WGfy 
W2Gfz=W2fz-WGfz 

* 

* CALCULATE VT2, VTG, AND VT2G 

* 


CALL 


CROSS (VT2fx,VT2fy,VT2fz,W2fx.W2fy,W2fz, 


Bfx.Bfy.Bfz) 


152 



CALL CROSS (VTGfx, VTGfy, VTGfz,WGfx,WGfy,WGfz,Bfx,Bfy,Bfz) 

VT2Gf x=VT2f x-VTGf x 
VT2Gf y=VT2f y-VTGf y 
VT2Gfz=VT2fz-VTGfz 

* 

* CALCULATE V(2G)GI AND V(2G)GII 

Vc 

CALL DOT (VGI , EGIf x , EGIfy , EGI f z , VT2Gf x , VT2Gf y , VT2Gf z) 

CALL DOT (VGII , EGI If x , EGI I f y , EGI If z, VT2Gf x , VT2Gf y , VT2Gf z) 

* 

* CALCULATE A13,A23,A33 

* 

CALL DET (DETI, W2Gf x , W2Gf y, W2Gf z, Nfx,Nfy,Nfz, EGI fx, EGIfy, EGI fz) 
A13=-KGI*VGI-DETI 

CALL DET (DETI I , W2Gf x , W2Gfy , W2Gf z, Nf x , Nf y , Nf z, EGIIf x , EGI I f y , EGI I f z) 
A23=-KGII*VGII-DETII 

CALL DET(DET3,Nfx,Nfy,Nfz,W2Gfx,W2Gfy,W2Gfz,VT2Gfx,VT2Gfy,VT2Gfz) 
CALL CROSS (Cx , Cy , Cz, W2f x , W2f y , W2f z , VTGf x , VTGfy , VTGfz) 

CALL CROSS (Dx , Dy , Dz, WGf x , WGf y , WGf z, VT2f x , VT2fy , VT2f z) 

CALL DOT (DET45 , Nf x , Nf y , Nf z , Cx~Dx , Cy-Dy , Cz~Dz) 
A33=KGI*VGI*VGI+KGII*VGII*VGII-DET3-DET45 

* 

* CALCULATE SIGMA 

P=A23*A23-A13*A13+ (KGI-KGI I) *A33 
SIGDBL=DATAN (2 . D00*A1 3*A23/P) 

SIGMA=0.5D00*SIGDBL 

* 

* CALCULATE K2I AND K2II 

* 

T1=P/ (A33*DCOS (SIGDBL) ) 

T2=KG I +KG 1 1 - ( A 1 3 *A 1 3 + A2 3 * A2 3 ) / A3 3 
K2I= (T1 + T2) /2.D00 
K2II=(T2-T1) /2.DOO 

Vc 

* CALCULATE E2I AND E2II 

Vc 

CALL ROTATE (E2lf x, E2I f y, E2I f z, EGIf x, EGIfy, EGI fz, -SIGMA, Nfx.Nfy, 

. Nfz) 

CALL ROTATE (E2I If x , E2I If y , E2I I f z, E2lf x , E2lf y , E2If z.RITAG, 

. Nfx,Nfy,Nfz) 

END 

* 

* FOR THE DETERMINATION OF PINION PRINCIPAL CURVATURES AND DIRECTIONS 

* 

SUBROUTINE PRINP1 (KS1 ,G1 , El If x , Ellfy , El If z, Elllf x , Elllfy , El Ilf z) 
IMPLICIT REAL*8(A-H,K,M-Z) 

REAL*8 TPAR (19) , LM 
CHARACTER *8 HG 
COMMON / T 1 / TP AR 
COMMON/AO/HG 

COMMON/A3/TND1 , TND2 , RITAG 

COMMON/B1/CSPH1 1 , SNPHl 1 , SP , EM, LM, CSRT1 , CSDl , SNDl , CSPSIP , SNPSIP 
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C0MM0N/B2/CSPIT1 , SNPIT1 ,MP1 ,MG2,QP 
C0MM0N/B3/Bf x , Bf y , Bf z 
C0MM0N/C2/Nf x , Nf y , Nf z 
COMMON/DI /UP , CSTAUP , SNTAUP 
KPI=CSPSIP/(UP*SNPSIP) 

KPII=O.DOO 
EPIf x=SNDl*SNTAUP 
EPIfy=CSTAUP 
EPIfz=CSDl*SNTAUP 

CALL ROTATE (EPIIf x , EPIIfy,EPIIfz,EPIfx, EPIfy , EPIf z, RITAG, 
. Nfx,Nfy, Nfz) 

* 

* CALCULATE W1P 

* 

IF (HG .EQ. ' L ' ) THEN 
Wlfx=-SNPIT1 
WPfx=-MPl*CSDl 
W1 fy=0. DOO 
WPfy=O.DOO 
Wlfz=-CSPIT1 
WPfz=MPl*SNDl 
ELSE 

Wlfx=SNPITl 
WPf x=MPl *CSD1 
W1 fy=0. DOO 
WPf y=0 . DOO 
Wlfz=CSPITl 
WPfz=~MPl*SNDl 
END IF 

WlPfx=Wlfx-WPfx 
WlPfy=Wl fy-WPfy 
WlPfz=Wlfz-WPfz 

* 

* CALCULATE VT2 , VTG, AND VT2G 

* 


CALL CROSS (VTlfx , VTl f y , VTl f z, W1 f x, Wlf y , Wlf z, Bf x , Bf y , Bfz) 
CALL CROSS (VTPlf x , VTP1 f y , VTP1 fz, WPf x , WPfy , WPfz, Bf x , Bfy , Bfz) 
IF (HG .EQ. 1 L 1 ) THEN 

CALL CROSS (VTP2fx,VTP2fy,VTP2fz,TPAR (18) ,EM,TPAR(19) , 

WPfx, WPfy, WPfz) 

ELSE 

CALL CROSS (VTP2fx,VTP2fy,VTP2fz,TPAR (18) , -EM, TPAR (19) , 
WPfx, WPfy, WPfz) 

END IF 

VTPf x=VTP 1 f x+VTP2f x 
VTPfy=VTPlfy+VTP2fy 
VTP f z=VTP 1 f z+VTP2f z 
VTlPfx=VTlfx-VTPfx 
VTlPf y= VTl fy- VTPf y 
VTlPf z = VTl fz~ VTPf z 

CALL DOT (VPI, EPIf x, EPIf y , EPI f z, VTIPfx , VTlPf y, VTlPf z) 

CALL DOT (VPII , EPI If x, EPIIf y.EPIIfz, VTIPfx, VTlPf y, VTlPf z) 
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* 

* CALCULATE A13,A23,A33 


CALL DET(DETI,WlPfx,WlPfy,WlPfz,Nfx,Nfy,Nfz,EPIfx,EPIfy,EPIfz) 

a13=-kpi*vpi-deti 

CALL DET (DETII ,WlPfx,WlPfy,WlPfz,Nfx,Nfy,Nfz, 
EPIIfx,EPIIfy,EPIIEz) 

A23=-KPII*VPII-DETII 

CALL DET (DET3 ,Nf x ,Nfy,Nfz,WlPfx,WlPfy,WlPfz, 

VTlPfx, VTlPfy , VTIPfz) 

CALL CROSS (Cx,Cy,Cz,Wlfx,O.DOO,Wlfz,VTPfx,VTy,z 

CALL CROSS (Dx ,Dy ,Dz, WPf x , 0 . D00 , WPf z, VTl f x, VTlfy , VTlf z) 

CALL DOT (DET45 ,Nf x ,Nfy , Nf z, Cx-Dx ,Cy Dy,Cz Dz; 

A33=KPI*VPI*VPI+KPII*VPII*VPII-DET3-DET45 


* 


* CALCULATE SIGMA 

P= A2 3 * A2 3 - A 1 3 * A 1 3+ (KP I -KP 1 1 ) * A3 3 
SIGDBL=DATAN (2 . D00*Al3*A23/P) 

SIGMA=0.5D00*SIGDBL 

* 

* CALCULATE KlI AND Kill 

* 

G1=P/ (A33*DCOS (SIGDBL) ) 

KS 1 =KP I +KP I I ~ (Al3*Al3 + A23*A23) / A33 
Kll= (KSl+Gl) /2.D00 
K1II=(KS1-G1)/2.D00 

* 

* CALCULATE Ell AND El I I 

CALL ROTATE (ElIfx,ElIfy, Ellfz, EPIfx,EPIfy,EPIfz, -SIGMA, Nfx.Nfy, 
CALL ROTATE(ElIIfx,ElIIfy,ElIIfz,ElIfx,ElIfy,ElIfz,RITAG, 

. Nfx,Nfy,Nfz) 

END 


I FOR THE DETERMINATION OF THE ANGLE BETWEEN GEAR PRINCIPAL DIRECTIONS 
* AND PINION PRINCIPAL DIRECTIONS 


* 


IBROUTINE SIGAN2(E2I£x,E2Ify,E2Ifz,E2IIfx,E2IIfy,E2lIfz,ElI x, 

•lIfy,ElIfz,CS2SIG,SN2SIG,SIGMl2) 

1PLICIT REAL*8(A-H,K,M-Z) 

VLL DOT(CSSIG,ElIfx,ElIfy,ElIfz,E2Ifx,E2lfy,E2Ifz) 

^LL DOT(SNSIG,ElIfx,ElIfy,ElIfz,-E2II£x, E2IIfy, E2II z) 

7GM2=4 . D00*DATAN (SNSIG/ (1 .DOO+CSSIG) ) 


SIGM12=.5D00*SIGM2 
CS2SIG=DCOS (SIGM2) 
SN2SIG=DSIN (SIGM2) 
END 


* 

* 

* 


FOR THE DETERMINATION OF CONTACT ELLIPS 

SUBROUTINE ELLIPS (KS2 , G2,KSl ,G1 , CS2SIG, SN2SIG.DEF , ALFAl , 
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AXISL, AXISS,ElIfx,ElIfy Ellfz) 

IMPLICIT REAL*8 (A-H.K.M-2) 

COMMON/ A3/TND1 , TND2,RITAG 

COMMON/C2/Nf x , Nf y , Nfz 

COMMON/E 1/XBf , YBf , ZBf 

D=DSQRT (G1*G1~2 . D00 ,V G1*G2 ,V CS2SIG+G2*G2) 

CS2AF1= (G1-G2’ V CS2SIG) /D 

SN2AF1=G2*SN2SIG/D 

ALFA1=DATAN (SN2AF1/ (1 .D00+CS2AFI) ) 

A= . 25D00 ,V DABS (KS 1-KS2-D) 

B= . 25D00*DABS (KS1-KS2+D) 

IF (KS2 . LT. KS1)THEN 
AXISL=DSQRT (DEF/A) 

AXISS=DSQRT (DEF/B) 

^LROTATCtXBt.VBE.ZBf.Enfx.Enfy.Enfz.R’WG-^j,^, 

ELSE 

AXISL=DSQRT (DEF/B) 

AXISS=DSQRT (DEF/A) 

CALL^ ROTATE (XBf.YBE.ZM.Enfx.SUfy.Ent^-^i.Nf^Nfy, 

END IF 

XBf =AXISL’ v XBf 
YBf =AXISL*YBf 
ZBf=AXISL*ZBf 
END 

COORDINATE TRANSFORMATION FOR F TO P 

SUBROUTINE PF (B2px , B2py , B2pz, Bf x , Bfy , Bfz) 

IMPLICIT REAL*8(A-H,K,M-Z) 

COMMON/ A4/ CSD2 , SND2 , CSPIT2 , SNP IT2 
C0MM0N/B4/CSPH2 , SNPH2 , CSPH21 , SNPH21 

[Mwf] = [Mwa] [Maf] 


CALL COMBI (wll,wl2,wl3,w21, w22 , w23 , w31 , w32, w33 wl w2 w3 

: S s S:6 s S:r°'- SN ™ 2, ’ CSPH2I ’ 0 - M0 ’ 0 --’° :D “--'“0, 

CALL TRCOOR (B2wx , B2wy , B2w Z , 

. wll , wl2, wl3, w21 , w22, w23, w31 , w32, w33,wl ,w2,w3, 

. Bfx.Bfy, Bfz) 

[Mpw] = [Mpa] [Maw] 

CALL COMBI (pll,pl2,pl3,p21, p22 , p23, P 31 , P 32, p 33 pi p2 0 3 

: o.M"o 0 Mri^o IT2 ’ 0 ' DOO,1 - DOO ' 0 - DOO - SNPIT2 ' 0 -“° :cs ” T2 ’ 

* O^DOoiofoGofo^DOO) '^^ PH2 C2PH2 ' °" D ^°’ 0 ' D ^ 0,0 " D ^^’ ^ 

CALL TRCOOR (B2px,B2py,B2pz, 
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. pll,pl2,pl3,p21,p22,p23,p31,p32,p33,pl ,p2,p3, 

. B2wx , B2wy , B2wz) 

END 

i< 

* USING EULER FORMULA TO DETERMINATION SURFACE INTERFERENCE 

* 

SUBROUTINE EULER (KS2 , G2 , KS1 , G1 , CS2SIG, SN2SIG, IEU) 

IMPLICIT REAL*8(A-H,K,M-Z) 

T=KS2-KS1 

U=DSQRT( (G2-G1*CS2SIG) **2+ (G1 ,V SN2SIG) **2) 

KR1=(T+U)/2.D00 
KR2=(T-U)/2.D00 
IF (KR1 *KR2 .LT. O.DOO)THEN 
IEU= 1 
ELSE 
IEU=0 
END IF 
END 

* 

* DETERMINANT 

i: 

SUBROUTINE DET (S , A, B , C , D, E , F , G, H, P) 

IMPLICIT REAL*8(A-H,K,M-Z) 

S=A*E*P+D*H*C+G*B*F-A*H*F-D*B*P-G*E*C 

RETURN 

END 

* 

* COORDINATE TRANSFORMATION 

* 

SUBROUTINE TRCOOR (XN , YN , ZN , R 1 1 , R 12 , R13 , R21 , R22 , R23 , R31 ,R32 , R33 , 
Tl,T2,T3,XP,YP,ZP) 

IMPLICIT REAL*8 (A-H,0-Z) 

XN=R11*XP+R12*YP+R13*ZP+T1 

YN=R21*XP+R22*YP+R23*ZP+T2 

ZN=R31*XP+R32*YP+R33*ZP+T3 

RETURN 

END 

* 

* MULTIPLICATION OF TWO TRANSFORMATION MATRICES 

Vf 

SUBROUTINE COMBI (Cl 1 , C12 , C13 , C21 , C22 , C23 , C31 , C32 , C33 , Cl , C2 , C3 , 
All ,A12,A13, A21 , A22.A23, A31,A32,A33,A1,A2,A3, 
B11,B12,B13,B21,B22,B23,B31,B32,B33,B1,B2,B3) 
IMPLICIT REAL*8(A-H,M,N,0-Z) 

C11=B31*A13+B21*A12+B11*A11 

C12=B32*A13+B22*A12+B12*A11 

C13=B33*A13+B23*A12+B13*A11 

C21=B31*A23+B21*A22+B11*A21 

C22=B32*A23+B22*A22+B12*A21 

C23=B33*A23+B23*A22+B13*A21 

C31=B31*A33+B21*A32+B11*A31 

C32=B32*A33+B22*A32+B12*A31 

C33=B33*A33+B23*A32+B13*A31 
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Cl=B3* * * * Al3+B2*Al2+Bl*All+Al 

C2=B3*A23+B2*A22+B1*A21+A2 

C3=B3*A33+B2*A32+B 1 *A3 1 +A3 

RETURN 

END 

* 

* DOT OF TWO VECTORS 

* 

SUBROUTINE DOT (V , XI , Y 1 , Z 1 , X2 , Y2 , Z2) 

IMPLICIT REAL*8(A-H,0-Z) 

V=X1*X2+Y1*Y2+Z1*Z2 

RETURN 

END 

* 

* CROSS OF TWO VECTORS 

* 

SUBROUTINE CROSS (X,Y,Z,A,B,C,D,E,F) 

IMPLICIT REAL*8(A-H,0-Z) 

X“B*F-C*E 

Y=C*D~A*F 

Z=A*E~B*D 

RETURN 

END 

* 

* ROTATION A VECTOR ABOUT ANOTHER VECTOR 

* 

SUBROUTINE ROTATE (XN , YN , ZN , XP , YP , ZP , THETA, UX , UY , UZ) 

IMPLICIT REAL*8(A-H,0-Z) 

CT=DCOS (THETA) 

ST=DSIN (THETA) 

VT=1 .DOO-CT 

Rll=UX ,v UX ,v VT+CT 

R12»UX*UY*VT-UZ*ST 

R13=UX*UZ*VT+UY*ST 

R 2 1 =UX*U Y * V T+UZ * S T 

R22-UY*UY*VT+CT 

R23=UY*UZ ,v VT-UX*ST 

R31=UX*UZ*VT~UY*ST 

R32-UY*UZ*VT+UX*ST 

R33=UZ*UZ*VT+CT 

CALL TRCOOR (XN, YN ,ZN,R11,R12,R13,R21 ,R22,R23 ,R31 ,R32, R33 , 

0 . D00 , 0 . D00 , 0 . D00 , 

XP , YP , ZP) 

RETURN 

END 

* 

* ***** SUBROUTINE NOLIN ***** 

* 

SUBROUTINE NONL IN (FUNC ,NSIG,NE,NC,X,Y,Y1, DELTA , A , I P VT , WORK) 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION X (NE) , Y (NE) , Y 1 (NE) , A (NE , NE) , IPVT (NE) , WORK (NE) 

EXTERNAL FUNC 

NDIM=NE 
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EPSI=1.D00/10.D00** * * NSIG 

CALL NONLIO (FUNC , EPS I , NE , NC , X , DELTA , NDIM , A, Y , Y 1 , WORK , I PVT) 

RETURN 

END 

* 

* ***** SUBROUTINE NOLINO ***** 

* 

SUBROUTINE NONLIO (FUNC , EPS I , NE , NC , X , DELTA, NDIM , A , Y , Y 1 , WORK , IP VT) 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION X(NE) ,Y(NE) ,Y1(NE) ,IPVT(NE) ,WORK(NE) ,A(NDIM,NE) 
EXTERNAL FUNC 

* NC: // OF COUNT TIMES 

DO 5 1=1, NC 
CALL FUNC (X, Y,NE) 

* NE: # OF EQUATIONS 

DO 15 J=1 , NE 

IF (DABS ( Y ( J) ) .GT.EPSI) GO TO 25 
15 CONTINUE 
GO TO 105 
25 DO 35 J=1 , NE 
35 Y1(J)=Y(J) 

DO 45 J= 1 , NE 

DIFF=DABS (X (j) ) *DELTA 

IF (DABS (X ( J) ) . LT. 1 . D-12) DIFF=DELTA 

XMAM=X ( J) 

X(J)=X(J)-DIFF 
CALL FUNC (X, Y, NE) 

X ( J) =XMAM 
DO 55 K= 1 , NE 

A(K, J) = (Y1 (K)-Y(K))/DIFF 
55 CONTINUE 
45 CONTINUE 
DO 65 J= 1 , NE 
65 Y(J)=-Y1(J) 

CALL DECOMP (NDIM,NE,A,COND, IPVT, WORK) 

CALL SOLVE (NDIM, NE, A, Y , IPVT) 

DO 75 J= 1 , NE 
X(J) =X(J)+Y(J) 

75 CONTINUE 
5 CONTINUE 
105 RETURN 
END 

3V 

* ***** SUBROUTINE DECOMP ***** 

* 

SUBROUTINE DECOMP (NDIM, N , A, COND, IPVT , WORK) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A (NDIM, N) , WORK (N) , IPVT (N) 

* 

* DECOMPOSES A REAL MATRIX BY GAUSSIAN ELIMINATION, 

* AND ESTIMATES THE CONDITION OF THE MATRIX. 
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* -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS - , BY G. E. FORSYTHE, 

* M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

Vc 

* USE SUBROUTINE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEM. 

Vc 

* INPUT.. 

* 

* NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A 

* N = ORDER OF THE MATRIX 

* A MATRIX TO BE TRIANGULAR I ZED 

* 

* OUTPUT . . 

& 

* A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PREMUTED 

* VERSION OF A LOWER TRIANGULAR MATRIX I - L SO THAT 

* (PERMUTATION MATRIX) *A=L*U 

* 

* COND = AN ESTIMATE OF THE CONDITION OF A. 

* FOR THE LINEAR SYSTEM A *X = B , CHANGES IN A AND B 

* MAY CAUSE CHANGES COND TIMES AS LARGE IN X. 

* IF COND+l.O .EQ. COND , A IS SINGULAR TO WORKING 

* PRECISION. COND IS SET TO 1.0D+32 IF EXACT 

* SINGULARITY IS DETECTED. 

* 

* IPVT = THE PIVOT VECTOR 

* IPVT (K) = THE INDEX OF THE K - TH PIVOT ROW 

* IPVT (N) = (-1)** (NUMBER OF INTERCHANGES) 

* 

* WORK SPACE. . THE VECTOR WORK MUST BE DECLARED AND INCLUDED 

* IN THE CALL. ITS INPUT CONTENTS ARE IGNORED. 

* ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. 

* 

* THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY 

* DET(A) = IPVT (N) * A(1 , 1) * A(2,2) * ... * A(N,N) . 

Vc 

IPVT (N) =1 

IF (N. EQ. 1) GO TO 150 
NM1=N - 1 

* COMPUTE THE 1-NORM OF A . 

AN0RM=0.D0 

DO 20 J=1,N 
T=0.D0 
DO 10 1=1, N 
10 T=T+DABS(A(I, J)) 

IF (T.GT.ANORM) ANORM=T 
20 CONTINUE 

* DO GAUSSIAN ELIMINATION WITH PARTIAL 

* PIVOTING. 

DO 70 K= 1 , NMl 

KP 1=K+ 1 

* FIND THE PIVOT. 

M=K 

DO 30 I=KPl ,N 


160 


IF (DABS (A (I , K) ) . GT . DABS (A (M,K) )) M=I 
30 CONTINUE 
IPVT(K) =M 

IF (M.NE.K) ipvt(n)=-ipvt(n) 

T=A (M, K) 

A(M,K)=A(K,K) 

A(K,K)=T 

SKIP THE ELIMINATION STEP IF PIVOT IS ZERO. 
IF (T.EQ.O.DO) GO TO 70 

COMPUTE THE MULTIPLIERS. 

DO AO I=KP 1 , N 
AO A(I,K)=-A(I,K)/T 

INTERCHANGE AND ELIMINATE BY COLUMNS. 

DO 60 J=KP1,N 
T=A (M, J) 

A(M, J) =a(k, J) 

A (K, J) =T ' 

IF (T.EQ.O.DO) GO TO 60 
DO 50 I=KP1 , N 

50 A (I , J) = A (I , J) +A (I , K) *T 

60 CONTINUE 
70 CONTINUE 


COND = (1-NORM OF A) * (AN ESTIMATE OF THE 1-NORM OF A-INVERSE) 

* THE ESTIMATE IS OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE 

* SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS 

* 0F EQUATIONS, (A-TRANSPOSE) *Y = E AND A*Z = Y WHERE E 

* IS A VECTOR OF +1 OR -1 COMPONENTS CHOSEN TO CAUSS GROWTH IN Y. 

* ESTIMATE = (1-NORM OF Z) / (1-NORM OF Y) 

it 

* SOLVE (A-TRANSPOSE) *Y = E . 

DO 100 K=1,N 

T=0 . DO 

IF (K.EQ.l) GO TO 90 
KM1=K-1 
DO 80 1=1, KM1 
80 T=T+A(I,K)*W0RK(I) 

90 EK= 1 . DO 

IF (T.LT.O.DO) EK=-1 . DO 

IF (A(K,K) .EQ.O.DO) GO TO 160 
Al 1=A(1 , 1) 

WORK (K) =- (EK+T) /A ( 1 , 1 ) 

100 CONTINUE 

DO 120 KB= 1 , NM1 
K=N-KB 
T=0. DO 
KP1=K+1 

DO 110 I=KP1,N 
110 T=T+A(l,K)*WORK(K) 
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WORK (K)=T 
M=IPVT(K) 

IF (M.EQ.K) GO TO 120 
T=WORK (M) 

WORK (M) “WORK (K) 

WORK (K) =T 

120 CONTINUE 

a 

YNORM=0 . DO 
DO 130 1=1, N 

130 YNORM=YNORM+DABS (WORK (I) ) 

* 

* SOLVE A*Z = Y 
CALL SOLVE (NDIM, N , A, WORK, IPVT) 

* 

ZNORM=0.D0 
DO 1 AO 1=1, N 

140 ZNORM=ZNORM+DABS (WORK (i) ) 

* 

* ESTIMATE THE CONDITION. 

COND= ANORM*ZNORM/ YNORM 

IF (COND.LT. 1 .DO) COND=1.DO 
RETURN 

* 1-BY-l CASE.. 

150 COND= 1 . DO 

IF (A(l, 1) .NE.O.DO) RETURN 

A 

* EXACT SINGULARITY 
160 C0ND=1 . 0D32 

RETURN 

END 

SUBROUTINE SOLVE (NDIM, N, A, B, IPVT) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION A(NDIM,N),B(N),IPVT(N) 

A 

* SOLVES A LINEAR SYSTEM, A*X = B 

* DO NOT SOLVE THE SYSTEM IF DECOMP HAS DETECTED SINGULARITY. 

Vc 

* -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS", BY G. E. FORSYTHE, 

* M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

* 

* INPUT.. 

A 

* NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A 

* N ORDER OF MATRIX 

a A TRIAN2ULARIZED MATRIX OBTAINED FROM SUBROUTINE DECOMP 

* B RIGHT HAND SIDE VECTOR 

A IPV T = PIVOT VECTOR OBTAINED FROM DECOMP 

A 

* OUTPUT.. 

A 
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B = SOLUTION VECTOR, X 

DO THE FORWARD ELIMINATION. 

IF (N.EQ.l) GO TO 50 
NM1=N-1 
DO 20 K=1 ,NM1 
KP1=K+1 
M=IPVT (K) 

T=B (M) 

B (M) =B (K) 

B(K)=T 

DO 10 I=KPl,N 
10 B(l)=B(l)+A(l,K)*T 

20 CONTINUE 

NOW DO THE BACK SUBSTITUTION. 

DO 40 KB=1 , NM1 
KM1=N-KB 
K=KM1+1 

B(K)=B(K) /A(K,K) 

T=-B(K) 

DO 30 1=1 , KM1 
30 B(I)=B(I)+A(I,K)*T 

40 CONTINUE 
50 B(1)=B(1)/A(1, 1) 

RETURN 

END 
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*********************************************************************** 

, * 

k 

* Gleason's Spiral Bevel Gears * 

* 

* Basic Machine-Tool Settings and Tooth Contact Analysis 

. k 

k 

* Curved Blade to Cut the Pinion * 

* * 
*********************************************************************** 

IMPLICIT REAL ,V 8(A-H,K,M-Z) 

REAL*8 X (1) , F (1) , FI (1) , PAR (6) ,LM,TX(6) ,TF(6) ,TF1 (6) ,TPAR(19) , 

AZSP (1,1) , WORKP (1) ,AZS(6,6) ,W0RK(6) , LANDAP , LANPDG , LANDPO 
CHARACTER*8 HG.HNGR 
DIMENSION IPVT(6) ,IPVTP(1) 

EXTERNAL PCN1 , PCN2, TCN 
COMMON/P 1 /PAR 
COMMON/T1/TPAR 
COMMON/AO/HG 

COMMON/Al/pl 1 ,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3 

COMMON/A2/SG,CSRT2,QG,SNPSIG,CSPSIG 

COMMON /A3/TND1 ,TND2,RITAG 

COMMON/A4/CSD2 , SND2 , CSPIT2 , SNPIT2 

COMMON/A5/CSQG , SNQG, THETAG 

COMMON/B 1 /CSPHl 1 , SNPH1 1 , SP , EM , LM , CSRT1 , CSD1 , SND1 , CSLANP , SNLANP 

COMMON/B2/CSPIT1 , SNPIT1 ,MP1 ,MG2 ,QP 

COMMON/ B 3/ B2fx , B2fy , B2fz 

C0MM0N/B4/ CSPH2 , SNPH2 , CSPH2 1 , SNPH21 

COMMON / B 5 / XCR , ZCR 

COMMON/C 1 /UG , CSTAUG , SNTAUG 

COMMON/C2/N2fx,N2fy,N2fz 

COMMON/D1/CSTAUP , SNTAUP 

COMMON /El/XBf,YBf,ZBf 

COMMON/ Fl/PHIGO 

COMMON/Gl/DAl ,DV1 

* 

* INPUT THE DESIGN DATA 


* 

k 

* TNI : number of pinion teeth 

* sec . 3.1 

* TN2 : number of gear teeth 

* sec . 3.1 

* RTldg , RTlmin : root angle of pinion (degree and arc minute, respec- 

* tively) 

* sec . 3.1 

* RT2dg, RT2min : root angle of gear (degree and arc minute, respec- 

* tively) 

* sec . 3.1 

* SHAFdg : shaft angle (degree) 

* sec. 3.1 

* BETAdg : mean spiral angle (degree) 

* sec. 3.1 

* ADIA : average gear cutter diameter 
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* sec. 3.1 

* W : point width of gear cutter 

* sec. 3.1 

* A : mean cone distance 

* sec. 3.1 

* ALPHdg : blade angle of gear cutter (degree) 

* sec. 3.1 

* RX : radius of blade 

* gear convex side 

* fig. 10 

* RV : radius of blade 

* gear concave side 

* fig. 10 

* DLTXdg : angle measured counterclockwise from root of gear to 

* the tangent of the contact path (degree) 

* gear convex side 

* fig. 19 

* DLTVdg : angle measured counterclockwise from root of gear to 

* the tangent of the contact path (degree) 

* gear concave side 

* fig. 19 

* M21XPR : first derivative of gear ratio 

* gear convex side 

* sec. 3.1.1 

* M21VPR : first derivative of gear ratio 

* gear concave side 

* sec. 3.1.1 

* AXILX : semimajor axis of contact ellipse 

* gear convex side 

* eq. (3.76) 

* AXILV : semimajor axis of contact ellipse 

* gear concave side 

* eq. (3.76) 

* HNGR : hand of gear ( ’ L 1 or * R 1 ) 

* DA : amount of shift along pinion axis 

* + : pinion mounting distance being increased 

* - : pinion mounting distance being decreased 

* DV : amount of pinion shaft offset 

* the same sense as yf shown in fig. 18 

* DEF : elastic approach 

* eq. (3.76) 

* EPS : amount to control calculation accuracy 

* 

* OUTPUT OF THE BASIC MACHINE-TOOL SETTINGS 

Vc 

* PSIGdg : gear blade angle 

* PSIPdg : pinion blade angle 

* RG : tip radius of gear cutter 

* RP : tip radius of pinion cutter 

* SG : gear radial 

* SP : pinion radial 

* QGdg : gear cradle angle 

* QPdg : pinion cradle angle 
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* MG2 : gear cutting ratio 

* MPl : pinion cutting ratio 

* EM : machining offset 

,c LM : machine center to back + sliding base 

,c XCR, ZCR : x and z coordinates of center of blade 

* 

DATA TNI , TN2/10.D00 ,41 .D00/ 

DATA RTldg,RTlmin/12.D00, 1.D00/ 

DATA RT2dg,RT2min/72.D00,25.D00/ 

DATA SHAFdg , BETAdg/90 . D00 , 35 . D00/ 

DATA ADIA/6.0D00/ 

DATA W/0.08D00/ 

DATA A/3.226D00/ 

DATA ALPHdg/ 20 . D00/ 

DATA DLTXdg/ 90.D00/ 

DATA DLTVdg/ 75.D00/ 

DATA M21XPR/-3.5D-03/ 

DATA M21VPR/5.2D-03/ 

DATA AXILX/0. 1710D00/ 

DATA AXILV/0. 1810D00/ 

DATA RX/40.00D00/ 

DATA RV/5O.O0D00/ 

DATA HNGR/'L'/ 

DATA DA,DV/0.D00, 0.D00/ 

DATA DEF/0. 00025DOO/ 

DATA EPS/l.D-12/ 

* 

* 

* 

DAl=DA 

DVl=DV 

HG=HNGR 

* 

* CONVERT DEGREES TO RADIANS 

* 

CNST=4 . D00*DATAN ( 1 . D00) / 1 80 . D00 

RITAG=90.D00 ,C CNST 

DLTX=DLTXdg*CNST 

DLTV=DLTVdg*CNST 

RT1= (RTldg+RTlmin/ 60. D00) *CNST 

RT2= (RT2dg+RT2min/ 60.D00) *CNST 

BETA=BETAdg*CNST 

PS IG=ALPHdg*CNST 

SHAFT=SHAFdg*CNST 

CSRT2=DCOS(RT2) 

SNRT2=DSIN(RT2) 

CSRTl=DCOS (RT1) 

SNRTl=DSIN (RT1) 

* 

* CALCULATE PITCH ANGLES 

* 

MM21=TN1/TN2 
c eq. (3. 1) 
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PITCH2=DATAN (DSIN (SHAFT) / (MM21+DC0S (SHAFT) ) ) 

IF(PITCH2 . LT. O.DOO) THEN 
PITCH2=PITCH2+180. DOO 
END IF 

CSPIT2=DC0S (PITCH2) 

SNPIT2=DSIN (PITCH2) 

c eq. (3.2) 

PITCH1=SHAFT-PITCH2 
CSPIT1=DC0S (PITCHl) 

SNPITl=DSIN (PITCHl) 

Vc 

* CALCULATE DEDENDUM ANGLES 

* 

c eq. (3.3) 

D1=PITCH1-RT1 

D2=PITCH2-RT2 

CSD1=DC0S(D1) 

SND1=DSIN(D1) 

TND1=SND1/CSD1 
CSD2=DCOS (D2) 

SND2=DSIN(D2) 

TND2=SND2/CSD2 

* 

* CALCULATE GEAR CUTTING RATIO 

* 

c e q . (3.7) 

MG2=DSIN (PITCH2) / CSD2 

* 

* FOR GEAR CONVEX SIDE 1=1, FOR GEAR CONCAVE SIDE 1=2. 

* 

DO 99999 1=1,2 
IF (I .EQ. 1)THEN 
WRITE(72,*) 'GEAR CONVEX SIDE' 

DLTA=DLTX 

M21PRM=M21XPR 

AXIL=AXILX 

R=RX 

ELSE 

WRITE (72,*) 'GEAR CONCAVE SIDE’ 

DLTA=DLTV 
M21PRM=M21VPR 
AXIL=AXILV 
R=RV 
END IF 
WRITE (72,*) 

c eq. (3.76) 

AXIA=DEF/ (AXIL* AXIL) 

* 

* CALCULATE GEAR BLADE ANGLE 

* 

c sec. 2.2 

IF (I .EQ. 2) THEN 
PSIG=180.D00*CNST-PSIG 
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END IF 

CSPSIG=DCOS (PSIG) 

SNPSIG=DSIN (PSIG) 

CTPSIG=CSPSIG/SNPSIG 

Vc 

* CALCULATE CUTTER TIP RADIUS 

Vr 

c eq. (3.8) 

IF (I .EQ. 1)THEN 
RG= (ADIA-W) /2.D00 
ELSE 

RG= (ADIA+W) / 2 . DOO 
END IF 

* 

* CALCULATE RADIAL 

* 

c eq. (3.9) 

IF (I .EQ. 1)THEN 

SG=DSQRT (ADIA*ADIA/4 . D00+A*A*CSD2*CSD2~A*ADIA*CSD2*DSIN (BETA) ) 

?’c 

* CALCULATE CRADLE ANGLE 

* 

c eq. (3.10) 

QG=DACOS ( (A*A*CSD2*CSD2+SG*SG-ADIA*ADIA/4.D00) / (2.D00*A*SG*CSD2) ) 
CSQG=DCOS (QG) 

SNQG=DSIN (QG) 

END IF 

* 

PAR (1)=RG*CTPSIG*CSPSIG 
PAR (4) =RG*CTPSIG 

* CALCULATE PHIG 

Vc 

PHIG=O.DOO 
PHIGO=PHIG 
CSPHIG=DCOS (PHIG) 

SNPHIG=DSIN (PHIG) 

5V 

IF (HG .EQ. ' L ' ) THEN 
IF (I .EQ. 1)THEN 

* Mmc=Mms*Msc 

c eq. (2.26) 

CALL COMB 1(11111,11112,11113, m2 1 , m2 2 , m2 3 ,m3 1 , m32 , m3 3 , ml , m2 , m3 , 
l.DOO,O.DOO,O.DOO,O.DOO,CSPHIG,SNPHIG,O.DOO,-SNPHIG,CSPHIG, 

0. DOO, 0. DOO, 0. DOO, 

1 . DOO , 0 . DOO , 0 . DOO , 0 . DOO , CSQG , -SNQG , 0 . DOO , SNQG , CSQG , 

0 . DOO , -SG*SNQG , SG*CSQG) 

END IF 

* Mpc=Mpm*Mmc 

c eqs. (2.25) , (3.13) 

CALL COMBI (pll,pl2,pl3,p21,p22,p23,p31, p32 ,p33,pl,p2,p3, 

CSD2, 0. DOO, -SND2 , O.DOO , 1 .DOO , O.DOO, SND2, 0.D00.CSD2 , 

0. DOO, 0. DOO, 0. DOO, 
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. ml 1 ,ml2 , ml3 ,m21 ,m22 ,m23,m31 ,m32 ,m33 ,ml ,m2,m3) 

* 

ELSE 

* 

IF (I .EQ. 1)THEN 

* Mmc=Mms*Msc 

c eq. (2.26) 

CALL COMBI (ml 1 ,ml2 ,ml3 ,m21 ,m22,m23,m31 ,m32,m33,ml ,m2,m3, 

1 . DOO , 0 . DOO , 0 . D00 , 0 . DOO , CSPHIG , -SNPHIG , 0 . DOO , SNPHIG, CSPHIG , 

0 . DOO , 0 . DOO , 0 . DOO , 

1 . DOO , 0 . DOO , 0 . DOO , 0 . DOO , CSQG , SNQG , 0 . DOO , -SNQG , CSQG , 

O.DOO , SG*SNQG, SG’''CSQG) 

END IF 

* Mpc=Mpm*Mmc 

c eqs. (2.25), (3.13) 

CALL COMBI (pll s pl2,pl3,p21 ,p22,p23,p31 ,p32,p33,pl,p2,p3, 

CSD2, 0. DOO, -SND2.0. DOO, 1. DOO, O.DOO, SND2, O.DOO, CSD2, 

O.DOO, O.DOO, O.DOO, 

. ml 1 ,ml2 ,ml3 , m2 1 ,m22,m23,m31 , m3 2 , m3 3 , ml , m2 ,m3) 

END IF 

* 

* DETERMINE MAIN CONTACT POINT 


* CALCULATE THETAG 

* 

c X ( 1 ) represents THETAG 

PAR (2) = (MG2-SNRT2) *CSPSIG 
IF (HG .EQ. ’ L * ) THEN 
PAR (3) =-SNQG*CSRT2*SNPSIG 

c step 1 in sec. 3.2 

X(1)=QG-BETA+RITAG 

ELSE 

PAR (3) =SNQG*CSRT2*SNPSIG 

c step 1 in sec. 3.2 

X(l)=- (QG-BETA+RITAG) 

END IF 

CALL NONLIN (PCNl , 14 , 1 , 100, X, F , FI , 1 . D-5 , AZSP , IPVTP , WORKP) 
THETAG=X(1) 

CSTHEG=DCOS (THETAG) 

SNTHEG=DSIN (THETAG) 

* 

* CALCULATE TAUG 

* 

c eq. (2.38) 

IF (HG .EQ. ' L 1 ) THEN 
TAUG=THETAG-QG+PHIG 
ELSE 

TAUG=THETAG+QG-PHIG 
END IF 

CSTAUG=DCOS (TAUG) 

SNTAUG=DSIN (TAUG) 
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* CALCULATE UG 

* 


c eq. (2.43) 

IF (HG .EQ. ' L ' ) THEN 

UG=RG*CTPS IG*CSPS IG-SG* ( (MG2-SNRT2) *CSPSIG*SNTHEG-DSIN (QG-PHIG) * 
// CSRT2*SNPSIG) / (CSRT2*SNTAUG) 

ELSE 

UG=RG*CTPSIG*CSPSIG-SG* ( (MG2-SNRT2) *CSPSIG*SNTHEG+DSIN (QG-PHIG) * 
// CSRT2*SNPSIG)/(CSRT2*SNTAUG) 

END IF 

* 

* CONVERT RADIAN TO DEGREE 

* 

PSIGDG=PSIG/CNST 
TAUGDG=TAUG/ CNST 
QGDG=QG/CNST 
THEGDG=THETAG/CNST 
PHIGDG=PHIGO/ CNST 

* 

* OUTPUT OF GEAR SETTINGS 

it 

WR I TE ( 7 2 , 10000) P S I GDG , QGDG , RG , SG , MG2 , TAUGDG , UG , THEGDG , PH I GDG 

* 

* CALCULATE MAIN CONTACT POINT 

* 

c eq. (2. 1) 

Bcx=RG*CTPSIG~UG*CSPSIG 

Bcy=UG*SNPSIG*SNTHEG 

Bcz~UG*SNPSIG*CSTHEG 

c eq. (2.2) 

Ncx=SNPSIG 

Ncy=CSPSIG*SNTHEG 

Ncz=CSPSIG*CSTHEG 

c eq. (2.9) 

EGIcx=0.D00 

EGIcy=CSTHEG 

EGIcz=“SNTHEG 

c eq. (3.13) 

CALL TRCOOR (Bpx , Bpy, Bpz, 

. pll ,pl2,pl3,p21,p22,p23,p31 ,p32,p33,pl,p2,p3, 

. Bcx,Bcy,Bcz) 

c eq. (3. 16) 

CALL TRCOOR (Npx , Npy , Npz, 

. pll,pl2,pl3»p21, p22 , p23 , p31 , p32 , p33, 0.D00, 0. D00, 0.D00 , 

. Ncx,Ncy,Ncz) 

c eq. (3.17) 

CALL TRCOOR (EGIpx , EGIpy, EGIpz, 

. pll,pl2,pl3,p21,p22,p23,p31,p32,p33,O.DOO,O.DOO,O.DOO, 

. EGIcx , EG Icy , EGIcz) 

c fig* 18 & sec. 3.3 

Bf x=Bpx 
Bf y=Bpy 
Bfz=Bpz 
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Nf x=Npx 
Nfy=Npy 
Nfz=Npz 
EGIf x=EGIpx 
EGIfy=EGIpy 
EGIfz=EGIpz 
>v 

* CALCULATE LANDAP 

* 

UGO=UG 

CSTAGO=CSTAUG 

SNTAGO=SNTAUG 

PHIGO=PHIG 

THETGO=THETAG 

* 

DO 99999 J=1 , 2 

* 

CSTAUG=CSTAGO 

SNTAUG=SNTAGO 

UG=UG0 

PHIG=PHIGO 

THETAG=THETGO 

IF ( J .EQ. 1)THEN 
WRITE (72,*) 'BLADE CONCAVE DOWN’ 

ELSE 

WRITE (72,*) 'BLADE CONCAVE UP' 

END IF 

LANDAP=DACOS(CSDl*Nfx-SNDl*Nfz) 

IF (I .EQ. 1)THEN 
IF (J .EQ. 1)THEN 
LANDAP=360 . DOO*CNST~ LANDAP 
PSIP=450.D00*CNST-LANDAP 
ELSE 

LANDAP= 180 . DOO*CNST-LANDAP 
PSIP=270.DOO*CNST- LANDAP 
END IF 
ELSE 

IF ( J .EQ. 1)THEN 
PS IP=90 . DOO*CNST-LANDAP 
ELSE 

LANDAP= 180 . D00*CNST+LANDAP 
PSIP=270 . D00*CNST- LANDAP 
END IF 
END IF 

CSLANP=DCOS (LANDAP) 

SNLANP=DS IN (LANDAP) 

* 

* CALCULATE TAUP 

Vc 

TAUP=DATAN2 (Nfy/SNLANP , (Nfx~CSDl*CSLANP) / (~SND1*SNLANP) ) 

IF ( J .EQ. 2) THEN 

TAUP=DATAN2 (-Nfy/SNLANP, (-Nf x~CSDl*CSLANP) / (-SNDl*SNLANP) ) 
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END IF 

CSTAUP=DCOS (TAUP) 

SNTAUP=DSIN (TAUP) 

* 

CALCULATE PRINCIPAL CURVATURES AND DIRECTIONS OF THE GEAR CUTTER 


c eq. (2.10) 

KGI=-CTPSIG/UG 
c eq. (2.12) 


KGII=0.D00 

c the second principal direction is determined by rotating of 

c the first principal derection about unit normal by 90 degrees 

CALL ROTATE(EGIIfx,EGIIfy,EGIIfz,EGIfx,EGIfy,EGIfz,RITAG 
. Nfx,Nfy,Nfz) 

* 

* CALCULATE W2G 

* 

c eqs. (3 . 18) - (3 . 20) 

IF (HG .EQ. ' L ' ) THEN 
W2fx=-SNPIT2 
WGfx=-MG2*CSD2 
W2fy=0 . D00 
WGf y=0 . D00 
W2fz=CSPIT2 
WGfz— MG2*SND2 
ELSE 

W2fx=SNPIT2 
WGf x=MG2*CSD2 
W2fy=0 . D00 
WGf y=0 . D00 
W2fz=-CSPIT2 
WGf z“MG2 *SND2 
END IF 

W2Gf x=W2f x-WGf x 
W2Gfy=W2f y-WGf y 
W2Gfz=W2fz~WGfz 

5'f 

* CALCULATE VT2 , VTG, AND VT2G 

* 

c eq. (3.22) 

CALL CROSS (VT2f x , VT2fy , VT2f z, W2f x , W2fy, W2fz, Bfx , Bf y , Bfz) 
c eq. (3.21) 

CALL CROSS (VTGf x , VTGf y , VTGf z, WGf x , WGfy , WGf z, Bf x , Bfy Bfz) 

c eq. (3.23) 

VT2Gf x=VT2f x-VTGf x 
VT2Gf y=VT2f y-VTGf y 
VT2Gfz=VT2f z-VTGfz 

Vc 

* CALCULATE V(2G)GI AND V(2G)GII 

Vc 

c eq . (3.24) 

CALL DOT (VGI , EGIf x , EGIfy , EGIfz, VT2Gfx, VT2Gfy, VT2Gfz) 
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c eq. (3.25) . 

CALL DOT (VGII , EGI If x , EGI I fy , EGI I f z, VT2Gf x , VT2Gfy , VT2Gf z) 

jV 

* CALCULATE A13,A23,A33 


- eq. (3.26) 

CALL DET(DETI,W2Gfx,W2Gfy,W2Gfz,Nfx,Nfy,Nfz,EGIfx,EGIfy,EGIfz) 
A13=-KGI*VGI-DETI 



A23=-KGII ,v VGII-DETII 
c eq. (3.28) 

CALL DET (DET3 , Nf x , Nf y , Nf z, W2Gf x , W2Gf y , W2Gf z, VT2Gf x , VT2Gf y , VT2Gf z) 
CALL CROSS (Cx , Cy ,Cz, W2f x , W2fy , W2fz, VTGfx, VTGfy, VTGfz) 

CALL CROSS (Dx , Dy , Dz , WGf x , WGf y , WGf z, VT2f x , VT2f y , VT2f z) 

CALL DOT (DET45 ,Nfx,Nfy,Nfz, Cx-Dx , Cy-Dy , Cz~Dz) 
A33=KGI*VGI*VGI+KGII*VGir v VGII-DET3-DET45 

* 

* CALCULATE SIGMA 

* 

c eq. (3.29) 

P«A23*A23-A13*A13+ (KGI-KGIl) *A33 
SIGDBL=DATAN (2 . D00*Al3*A23/P) 

SIGMA=0.5D00*SIGDBL 

* 

* CALCULATE K2I AND K2II 

* 

c eqs . (3 . 30) - (3 . 31) 

T1=P/ (a33*DCOS (SIGDBL) ) 

T2=KGI+KGI I - (Al3 ,: Al3+A23*A23) / A33 
K2I= (T1+T2) /2.D00 
K2II=(T2-T1)/2.D00 

:‘c 

* CALCULATE E2I AND E2II 

Vc 

c description after eq. (3.29) 

CALL ROTATE (E2lf x,E2lfy ,E2If z, EGIf x , EGIfy , EGIfz, “SIGMA, Nfx.Nfy, 

. Nfz) 

CALL ROTATE (E2IIfx,E2l If y,E2Hfz,E2Ifx,E2lfy,E2Ifz,RITAG, 

. Nfx.Nfy, Nfz) 
c eq. (3.44) 

TNETAG=DSIN (DLTA+SIGMA) /DCOS (DLTA+ SIGMA) 

* 

* CALCULATE W2 

* 

c eq. (3.33) 

IF (HG .EQ. 1 L ' ) THEN 
W2fx=-MM21*SNPIT2 
W2fy=0.D00 
W2fz=MM21*CSPIT2 
ELSE 

W2fx=MM21*SNPIT2 

W2fy=0.D00 
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W2fz=-MM21*CSPIT2 
END IF 


* CALCULATE W1 

* 

c eq. (3.32) 

IF (HG .EQ. ' L ' ) THEN 
Wlfx=-SNPIT1 
Wlfy=O.DOO 
Wlfz=~CSPITl 
ELSE 

Wlf x=SNPITl 
Wlfy-O.DOO 
Wlfz=CSPITl 
END IF 

* 

* CALCULATE W12 

* 

c eq. (3 . 34) 

W12fx=Wlfx-W2fx 

W12fy=Wlfy-W2fy 

W12fz=Wlfz-W2fz 

* 

* CALCULATE VT2 

* 

c eq. (3.36) 

CALL CROSS (VT2fx , VT2fy, VT2fz, W2f x, W2fy, W2fz, Bf x,Bfy , Bfz) 

vt 

* CALCULATE VT1 

* 

c eq. (3.35) 

CALL CROSS (VTlfx, VTlfy, VTlfz, Wlfx, Wlfy,Wlfz,Bfx,Bfy,Bfz) 

* CALCULATE VT12 

* 

c eq. (3.37) 

VT12fx=VTlfx-VT2fx 

VT12fy=VTlfy-VT2fy 

VT12fz=VTlfz-VT2fz 

* 

* CALCULATE V2 

* 

c eq. (3.38) 

CALL DOT (V2l,VT12fx, VT12fy, VT12fz, E2lfx, E2lfy, E2lfz) 
c eq. (3.39) 

CALL DOT (V2lI,VT12fx,VT12fy,VT12fz,E2lIfx, E2llf y , E2I If z) 

* 

* CALCULATE A31 

* 

c eq. (3.40) 

CALL DET (DET1 ,W12fx,Wl2fy,W12fz,Nfx,Nfy,Nfz,E2lfx,E2lfy,E2lfz) 
A31=-K2I*V2I-DET1 
c eq. (A. 33) 
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A13=A31 


* CALCULATE A3 2 

* 

c eq. (3.41) 

CALL DET(DET2,Wl2fx,W12fy,W12fz,Nfx,Nfy,Nfz,E2IIfx,E2IIfy,E2IIfz) 
A32=-K2II*V2II-DET2 

c eq. (A. 35) 

A23=A32 

* 

* CALCULATE A3 3 

* 

c eq. (3 . 42) 

CALL DET (DET3,Nfx,Nfy,Nfz, W12f x, W12fy,Wl2fz, VT12fx, VT12fy, VTl2fz) 
CALL CROSS (Cx, Cy , Cz, Wlf x , W1 fy , Wlf z, VT2f x , VT2fy , VT2fz) 

CALL CROSS (Dx , Dy ,Dz, W2f x , W2fy , W2f z, VTlf x , VTl fy, VTlfz) 

CALL DOT (DOT l,Nfx,Nfy,Nfz, Cx-Dx , Cy~Dy , Cz _ Dz) 

CALL DET (DET4 , Nf x , Nf y , Nf z, W2f x ,W2fy,W2fz,Bfx,Bfy,Bfz) 
A33=K2I*V2I*V2I+K2II*V2 II’ v V2II-DET3-DOT1+M21PRM*DET4 

* 

* CALCULATE ETAP 

* 

c eq. (3.53) 

ETAP=DATAN ( ( (A33+A31*V2I) *TNETAG-A31*V2ll) / (A33+A32* 

. (V2II-V2I*TNETAG))) 

TNETAP=DSIN (ETAP) /DCOS (ETAP) 

* 

* CALCULATE All, Al2, AND A22 

* 

N3= (1 . DOO+TNETAP*TNETAP) *A33 

c eq. (3.72) 

N1=(A13*A13-(A23*TNETAP)**2)/N3 
c eq. (3.73) 

N2= (A23+A1 3*TNETAP) * (Al3+A23*TNETAP) /N3 

KS2=K2I+K2II 

G2=K2I-K2II 

c eqs . (3.74) , (3.75) 

KS1=KS2- ( (4.D00*AXIA*AXIA-N1*N1-N2*N2) * ( 1 . DOO+TNETAP*TNETAP) / 

. (-2 . D00*AXIA* (1 . DOO+TNETAP*TNETAP) +N1* (TNETAP*TNETAP~1 . D00) 

. -2.D00*N2*TNETAP)) 

c eqs. (3.66), (3.69) & description after eq. (3.60) 

A11=TNETAP*TNETAP/ (1 . DOO+TNETAP*TNETAP) * (KS2-KS1) -*-N 1 

c eqs. (3.67), (3.70) & description after eq. (3.60) 

A12=-TNETAP/ (1 . DOO+TNETAP*TNETAP) * (KS2-KS1) +N2 


c eqs. (3.68), (3.71) & description after eq. (3.60) 

A22= 1 . D00 / (1 . DOO+TNETAP*TNETAP) * (KS2-KS 1) ~N1 

c eq. (A. 32) 

A21=A12 


* CALCULATE SIGMA (12) 

* 

c eq. (3.77) 

SIGDBL=DATAN (2 . D00*Al2/ (K2I-K2II-A1 1+A22) ) 
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SIGM12=.5D00*SIGDBL 


* CALCULATE K1I AND Kill 

* 

c eq. (3.78) 

G1=2.D00*A12/DSIN (SIGDBL) 

c eq. (3.79) 

Kl I=. 5D00* (KS1+G1) 

Kl 11= . 5D00* (KS1-G1) 

* 

* CALCULATE Ell AND El I I 

Vc 

c similar to description after eq. (3.29) 

CALL ROTATE (ElIfx,ElIfy,ElIfz,E2lfx,E2lfy,E2lfz, - SIGM12,Nfx,Nfy, 
. Nfz) 

CALL ROTATE (ElIIfx,ElIIfy,ElIIfz,ElIfx,ElIfy,ElI fz, RITAG, 

. Nfx,Nfy,Nfz) 

* 

* PINION 


* CALCULATE PRINCIPAL DIRECTIONS OF THE PINION CUTTER 

* 

c eq. (3.92) 

IF (HG .EQ. ' L ' ) THEN 
EPIfx=SNDl*SNTAUP 
EPIfy=CSTAUP 
EPIf z=CSDl*SNTAUP 
ELSE 

EPIfx=-SNDl*SNTAUP 
EPIfy=-CSTAUP 
EP I f z=-CSDl *SNTAUP 
END IF 

IF (DACOS (EGIfx*EPIfx+EGIfy*EPIfy+EGIfz*EPIfz)/CNST .GT. 45.D00) 

. THEN 

EPIfx=-EPIfx 
EPIfy=-EPIfy 
EPIfz=-EPIfz 
END IF 

* 

CALL ROTATE (EPI If x, EPI I fy, EPI I fz, EPIf x, EPIf y, EPIf z, RITAG, 

. Nfx,Nfy,Nfz) 

* 

* CALCULATE THE ANGLE BETWEEN PRINCIPAL DIRECTIONS OF PINION AND CUTTER 


c cross product of eli and epi 

SNSIGM=(ElIfy*EPIfz-ElIfz*EPIfy)/Nfx 
c dot product of eli and epi 


CSSIGM=E1 If x*EPIf x+EHfy*EPIfy+ElIf z*EPIfz 
CS2SIG=2.D00 ,V CSSIGM*CSSIGM-1.D00 
TN2SIG=2 . D00*SNSIGM*CSSIGM/CS2SIG 

* 

* CALCULATE PRINCIPAL CURVATURES OF PINION CUTTER 
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c eq. (2.20) 

KPII=1 .D00/R 
IF ( J .EQ. 2) THEN 
KPII=-KPII 
END IF 

c eq. (3.94) 

KPI=(KPir v (KlI*CSSIGM*CSSIGM+KlII*SNSIGM*SNSIGM)-KlI w KlIl)/ 

. (KPII-Kl I*SNSIGM*SNSIGM-Kl I I*CSSIGM Vt CSSIGM) 

Vc 

’’'CALCULATE All, Al2, AND A22 
c eq. (A. 31) 

All=KPI-KlI’''CSSIGM’' ! CSSIGM-KlII’''SNSIGM*SNSIGM 

c eq. (A. 32) 

A12= (Kll-Klll) *SNSIGM*CSSIGM 
c eq. (A. 34) 

A22=KPII-Kir v SNSIGM’ , 'SNSIGM-KlII*CSSIGM’ , 'CSSIGM 

* 

* CALCULATE ZCR 

* 

c eq. (3.101) 

IF ( J .EQ. 1)THEN 
ZCR= (SNLANP/KPI) -R*SNLANP 
ELSE 

ZCR=- (SNLANP/KPI) -R*SNLANP 
END IF 

* 

* CALCULATE XCR 

* 

c eq. (3.99) 

Bmx=-Bfx’ v CSDl+Bfz’''SNDl 

XCR=Bmx-R’ c CSLANP 

* CALCULATE RP 

c eq. (3.103) 

IF (I* J .EQ. 2) THEN 
RP=ZCR+DSQRT (DABS (R*R-XCR’' C XCR) ) 

ELSE 

RP=ZCR-DSQRT (DABS (R’ V R-XCR’ , 'XCR) ) 

END IF 

* 

* CALCULATE MCP 

"Sc 

Zll=Nfy*EPIfz-Nfz*EPIfy 

Z12=Nfy*EPIfx-Nfx’ , 'EPIfy 

Z21=Nfy*EPIIfz-Nfz*EPIIfy 

Z22=Nfy*EPIIfx-Nfx*EPIIfy 

c eqs . (3.107), (3.108) 

Cl l=Zl 1*CSD1+Z12*SND1 

C12=-Z11*SNPIT1+Z12*CSPIT1 

C22=-Z21*SNPIT1+Z22*CSPIT1 
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IF (HG .EQ. ' R ' ) THEN 
Cl 1=-C1 1 
C12—C12 
C22=-C22 
END IF 

c eq. (3.119) 

T4= (Bfy*CSRTl) / (EPIIfx*CSDl-EPIIf z *SNDl) 

IF (HG .EQ. ' R ' ) THEN 
T4=-T4 
END IF 

c eq. (3.120) 

T1=-C11/KPI 

T2= (A1 1 *KP I I *T4+A1 1 *C22-A12*C 1 2) / (Al2*KPl) 

c eq. (3.122) 

Ull-Tl*EPIfx 
U12=T2*EPIf x+T4*EPIIfx 
U21=Tl*EPIfy 
U22=T2*EPIfy+T4*EPIIfy 
U31=Tl*EPIfz 
U32»T2*EPIfz+T4*EPIIfz 
c eq. (3.124) 

V1=U21* (Nfz*CSDl+Nfx*SNDl) -Nf y* (Ul 1*SND1+U31*CSD1) 
c eq. (3.125) 

V2= (U22 X CSD1-U21*SNPIT1) *Nfz~ (Ul l ,t CSPITl+U12*SNDl+U32 A CSDl-U31 
*SNPIT1) 5 " t Nfy+ (U21*CSPITl+U22*SNDl)*Nfx 
c eq. (3.126) 

V3=U22*CSPITl*Nfx+(U32*SNPITl-U12*CSPITl)*Nfy-U22 A SNPITl*Nfz 

IF (HG .EQ. 1 R ' ) THEN 
V1=-V1 
V 2=-V2 
V3=-V3 
END IF 

c eq. (3.132) 

HI 1 = -U21*CSPIT1+SND1* (Bfz ,v SNPITl-Bfx*CSPITl) 

c eq. (3.134) 

H21=Ull*CSPITl-U31*SNPITl+Bfy*SNRTl 
c eq. (3.136) 

H31=U21*SNP IT1+CSD1* (Bfz*SNPITl-Bf x*CSPITl) 
c eq. (3.133) 

H12= (Bfz Vr SNPITl-Bfx*CSPITl-U22) *CSPIT1 

c eq. (3.135) 

H22=-(Bfy-U12*CSPITl+U32*SNPITl) 
c eq. (3.137) 

H32=- (Bfz*SNPITl-Bfx*CSPITl-U22) *SNPIT1 
IF (HG .EQ. 1 R ' ) THEN 

Hll=U21*CSPITl+SNDl*(Bfz*SNPITl-Bfx*CSPITl) 

H21— Ull*CSPITl+U31*SNPITl+Bfy*SNRTl 
H31=-U21*SNPITl+CSDl*(Bfz*SNPITl-Bfx*CSPITl) 

H12= (Bfz*SNPITl-Bfx*CSPITl+U22) *CSPIT1 
H22=-(Bfy+U12*CSPITl-U32*SNPITl) 

H32=- (Bfz*SNPITl-Bfx*CSPITl+U22) *SNPIT 1 
END IF 

c eq. (3.139) 
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Fl-Nfx*Hll+Nfy*H21+Nfz*H31 

c eq. (3.140) 

F2=Nfx*H12+Nfy*H22+Nfz*H32 

c eq. (3.145) 

Y2=A12*(2.D00*KPI*T1*T2-V2-F1) 

wi-Y3/Y2 I * T2 * T2tKPIP>I '‘* T4 ' V3 ' F2) ’ <KPI * T2+C12)< ' <KPII * T4tC22) 

* 

* CALCULATE EM AND LM 

* 

c eq. (3.122) 

VTlPfx=Ull*MPl+U12 
VT 1 P f y =u2 1 *MP 1 +U22 
VTlPf Z =U31*MP1+U32 

c eq. (3.111) 

IF(HG .EQ. ' L ' ) THEN 

EM= (Bfy*CSPITl-VTlPfx)/(MPl*SNDl)+Bfy 

LM- (B fx* cspiT1 - Bfz * sNpiTi +VT ip fy) / MP i +Bfx *SNDH-Bfz*CSDl 
ELSE 

EM= (-Bfy*CSPITl-VTlPfx) / (MP1*SND1)-Bfy 
^ = ^ fX * CSPIT1 " Bfz ’'' SNPin_VT1Pf y ) / MP1+Bfx ’ V SNDl + Bfz’ f 'CSDl 

END IF 

* 

* CALCULATE SP AND QP 

* 

c eqs. ( 3-150 ) _ ( 3>151) 

IF(HG .EQ. ’ L ’ ) THEN 
IF ( J .EQ. 1)THEN 
Zl=-Bfy+EM-SNLANP/KPI ,V SNTAUP 
Z2=Bfx*SNDl+Bfz*CSDl-LM-SNLANP/KPI*CSTAUP 

ELSE 

Zl=-Bfy+EM+SNLANP/KPI ,V SNTAUP 

Z2=Bfx*SNDl+Bfz*CSDl-LM+SNLANP/KPI*CSTAUP 

END IF 
ELSE 

IF ( J .EQ. 1)THEN 
Z 1 =B f y+EM+SNLANP /KP I *SNTAUP 

Z2=Bfx*SNDl+Bfz*CSDl-LM-SNLANP/KPI*CSTAUP 

ELSE 

Z 1 =B f y+EM-SNLANP/KP I *SNTAUP 

Z2=Bfx*SNDl+Bfz*CSDl-LM+SNLANP/KPI*CSTAUP 

END IF 
END IF 

SP=DSQRT (Z1*Z1+Z2’ V Z2) 

QP=DATAN(Z1/Z2) 

IF (HG .EQ. ' L 1 ) THEN 

thetap=taup-qp 

ELSE 

THETAP=TAUP+QP 
END IF 

* 

* CONVERT RADIAN TO DEGREE 

* 
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PSIPDG=PSIP/CNST 
TAUPDG=TAUP/CNST 
QPDG=QP/CNST 
THEPDG= THET AP / CN S T 
LAN P DG= LANDAP/CNST 


* OUTPUT 


WRITE (72, 10001) PSIPDG,QPDG, 

THEPDG 

10000 FORMAT (IX, ’GEAR SETTINGS:',/ 

, IX , ' PSIGDG = ' ,G20. 12, 12X, 

, IX, ' RG =' ,G20. 12, 12X, 

, IX, ' MG2 =' ,G20. 12, 12X, 

, IX, ' UG = ' ,G20. 12, 12X, 

’ IX, 'PHIGODG = ' ,G20. 12, // 

, IX, 'PINION SETTINGS:',/) 
10001 ’ FORMAT (IX, 1 PSIPDG = ’ ,G20. 12, 12X, 

IX »pp = ,G20. 12, 1/A, 

IX ’ MPl = 1 , G20 . 12 , 1 2X , 

’.IX/XCR = ' , G20 . 12 , 12X, 

, IX, 'EM =' ,G20. 12, 12X, 

, IX , ' TAUPDG =' ,G20. 12, 12X, 


RP , SP ,MP 1 , LANPDG.XCR , ZCR , EM, LM, TAUPDG, 


' QGDG 
’SG 

' TAUGDG 
' THETAGDG 


'QPDG 

'SP 

' LANDAPDG 

'ZCR 

'LM 

' THETAPDG 


>' , G20. 12, / 
■' ,G20. 12, / 
■' .G20.12,/ 
=' ,G20. 12,/ 


= ’ ,G20. 12,/ 

=' ,G20. 12,/ 
=’ ,G20. 12, / 
■' ,G20. 12,/ 
=' , G20. 12, / 

= ' ,G20. 12, /) 


* TCA 

TPAR (1) -RG*CSPSIG/SNPSIG*CSPSIG 
TPAR (2) = (MG2-SNRT2) *CSPSIG 
TPAR (3) =CSRT2*SNPSIG 
TPAR (4) =RG ,V CSPSIG/ SNPSIG 

TPAR(5)=CSD2*SNPSIG 
TPAR (6) =SND2 ,c CSPSIG 
TPAR (7) =SND2*SNPSIG 
TPAR (8) =CSD2*CSPSIG 
TPAR(9 )=ZCR ,v CSRT 1 
TPAR (10) =SP*CSRT1 
TPAR (11) =EM*CSRTl 
TPAR (12) =XCR*CSRTl 
TPAR (13) =SP ,V (MP1-SNRT1) 
TPAR(14)=EM*SNRT1 
TPAR(15)=LM*SNRT1 
TPAR (16) =J 
TPAR (17) =R 

* 

PHIP=0 . D00 
PHI21=0.D00 
PHI11=0.D00 
CSPH1 l=DCOS (PHIl 1) 

SNPH11 = DSIN (PHIll) 


TX(1)=PHIP 
TX (2) =THETAP 
TX(3) =LANDAP 


180 



TX(4)=PHI21 
TX (5) =PHIG 
TX(6)=THETAG 

CALL NONLIN (TCN , 1 4 , 6 , 200 , TX , TF , TF 1 , 1 . D~5 , AZS , IPVT , WORK) 
PHIP0=TX(1) 

THEP0=TX(2) 

LANDP0=TX (3) 

PHI210=TX (4) 

PHIG0=TX (5) 

THEG0=TX (6) 

TX(1) =PHIP0 
TX(2)=THEP0 
TX (3) =LANDP0 
TX (4) =PHI210 
TX(5)=PHIG0 
TX (6) =THEG0 

DPHI 1 1= 18 . D00/36 . D00*CNST 

t 

DO 100 I J= 1 , 60 
CSPHl l=DCOS (PHI 11) 

SNPH1 1=DSIN (PHI 1 1) 

CALL NONLIN (TCN , 1 4 , 6 , 200 , TX , TF , TF 1 , 1 . D-5 , AZS , IPVT , WORK) 

PHIP=TX (1) 

THETAP=TX (2) 

LANDAP=TX (3) 

PHI21=TX (4) 

PHIG=TX (5) 

THETAG=TX (6) 

ERR OR = ( (PHI 2 1 *36 . D02-PH 1 2 1 0*36 . D02) -PHI 1 1 *36 . D02*TN 1/TN2) /CNST 

J*C 

CALL PRING2(KS2,G2,E2lfx,E2Ify,E2lfz,E2lIfx,E2lIfy,E2IIfz) 

CALL PRINPl (KSl,Gl,ElIfx,ElIfy,ElIfz,ElIIfx,ElIIfy,ElIIfz) 

CALL SIGAN2(E2lfx,E2Ify,E2Ifz,E2IIfx,E2IIfy,E2IIfz,ElIfx,ElIfy, 
ElI£z,CS2SIG,SN2SIG,SIGMl2) 

CALL EULER (KS2 , G2 , KS 1 , G1 , CS2SIG, SN2SIG, IEU) 

IF (IEU .EQ. 1)THEN 
WRITE (72,*) ’THERE IS INTERFERENCE’ 

GO TO 88888 
END IF 

* 

CALL ELL I PS (KS2 ,G2,KS1,G1,CS2SIG,SN2SIG, DEF , ALFAl , 

AXISL, AXISS,ElIfx,ElIfy,ElIfz) 

Vf 

CALL PF(B2px,B2py,B2pz,B2fx,B2fy,B2fz) 

* 

* XBf , YBf , and ZBf is the direction of the long axis of the ellipse 

* 

CALL PF (XBp ,YBp,ZBp,XBf,YBf,ZBf) 

ELBlpx=B2px+XBp 

ELBlpz=B2pz+ZBp 

ELB2px=B2px-XBp 

ELB2pz=B2pz-ZBp 
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IF(I .EQ. 1 .AND. J .EQ. 1 ) THEN 
WRITE (9,9000) IJ,PHI11/CNST, IJ, ERROR 
IF(IJ .LE. 37) THEN 
WRITE (8 , 8000) I J , B2pz, IJ,B2px 
WR ITE (7 , 7000) ELB lpz , ELB lpx , ELB2pz, ELB2px 
END IF 

ELSE IF (I .EQ. 1 .AND. J .EQ. 2) THEN 
WRITE (79, 9000) IJ,PHI11/CNST, I J, ERROR 
IF (I J .LE. 37) THEN 
WRITE (78 , 8000) IJ,B2pz, I J , B2px 
WRITE (77 , 7000) ELBlpz, ELBlpx, ELB2pz,ELB2px 
END IF 

ELSE IF (I .EQ. 2 .AND. J .EQ. 1 ) THEN 
WRITE (29 , 9000) I J , PHI 1 1 /CNST , I J , ERROR 
IF (I J .LE. 37) THEN 
WRITE (28, 8000) IJ,B2pz, IJ,B2px 
WRITE (27 , 7000) ELBlpz, ELBlpx , ELB2pz, ELB2px 
END IF 
ELSE 

WRITE (89 , 9000) I J , PHI 1 1 /CNST , I J , ERROR 
IF (I J .LE. 37) THEN 
WRITE (88 , 8000) IJ,B2pz, I J , B2px 
WRITE (87 , 7000) ELBlpz, ELBlpx , ELB2pz,ELB2px 
END IF 
END IF 

* 

PH II 1=PHI 1 1+DPHI 1 1 

* 

100 CONTINUE 


PHI 1 1=0. D00 

TX(1)=PHIP0 
TX (2) =THEP0 
TX (3) =LANDP0 
TX (4) =PHI210 
TX (5) =PHIG0 
TX (6) =THEG0 

DO 200 I J=1 ,60 
CSPH1 l=DCOS (PHI1 1) 

SNPH1 1=DSIN (PHI 1 1) 

CALL NONLIN (TCN, 1 A , 6 , 200 , TX, TF , TF1 , 1 .D~5, AZS, IP VT, WORK) 

PHIP=TX (1) 

THETAP=TX (2) 

LANDAP=TX (3) 

PHI21=TX (A) 

PHIG=TX (5) 

THETAG=TX (6) 

ERROR= ( (PHI21 *36 . D02-PHI2 10*36 . D02) -PHI 11*36. D02*TN1/TN2) /CNST 
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CALL PRING2 (KS2 ,G2,E2Ifx,E2lfy,E2Ifz,E2lIfx,E2IIfy,E2IIfz) 

CALL PRINPl (KS1 ,Gl,ElI£x,ElIfy,ElIfz,ElIIfx,ElIIfy,ElIIfz) 

CALL SIGAN2(E2Ifx,E2Ify,E2Ifz,E2lIfx,E2lIfy,E2IIfz,ElIfx,ElIfy, 
ElIfz,CS2SIG,SN2SIG,SIGMl2) 

CALL EULER (KS2 , G2 , KS 1 , G1 , CS2S IG , SN2SIG, IEU) 

• IF (IEU .EQ. 1)THEN 

WRITE (72,*) 'THERE IS INTERFERENCE' 

GO TO 88888 
END IF 

* 

CALL ELLIPS (KS2 , G2 , KS1 , G1 , CS2SIG, SN2SIG, DEF , ALFAl , 

AXISL , AXISS ,ElIfx,ElIfy,ElIfz) 

* 

CALL PF(B2px,B2py,B2pz,B2fx,B2fy,B2fz) 

* 

* XBf, YBf, and ZBf is the direction of the long axis of the ellipse 

Vr 

CALL PF (XBp , YBp , ZBp , XB£ , YBf , ZBf ) 

ELBlpx=B2px+XBp 

ELBlpz=B2pz+ZBp 

ELB2px=B2px-XBp 

ELB2pz=B2pz-ZBp 

* 

IF (I .EQ. 1 .AND. J .EQ. 1) THEN 
WRITE(9, 9001) U.PHIll/CNST, IJ, ERROR 
IF (I J .LE. 37) THEN 
WRITE (8 , 8001) IJ,B2pz,IJ,B2px 
WRITE (7 , 7000) ELBlpz, ELBlpx , ELB2pz, ELB2px 
END IF 

ELSE IF (I .EQ. 1 .AND. J .EQ. 2) THEN 
WRITE(79,9001)IJ,PHIll/CNST, I J, ERROR 
IF (I J .LE. 37) THEN 
WRITE (78 , 8001) I J , B2pz, I J , B2px 
WRITE (77 , 7000) ELBlpz, ELBlpx , ELB2pz, ELB2px 
END IF 

ELSE IF (I .EQ. 2 .AND. J .EQ. 1) THEN 
WRITE (29 , 9001) I J , PHI1 1/CNST , I J, ERROR 
IF (I J .LE. 37) THEN 
WRITE (28, 8001) IJ,B2pz, I J,B2px 
WRITE (27, 7000) ELBlpz, ELBlpx,ELB2pz,ELB2px 
END IF 
ELSE 

WRITE(89, 9001) IJ.PHIll/CNST.IJ, ERROR 
IF (I J .LE. 37) THEN 
WRITE (88, 8001) IJ,B2pz, IJ,B2px 
WRITE (87 , 7000) ELBlpz, ELBlpx , ELB2pz, ELB2px 
END IF 
END IF 

* 

PHIl 1=PHI 1 1-DPHI 1 1 
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CONTINUE 



99999 CONTINUE 

88888 CONTINUE 

?v 

7000 FORMAT (6X, ' EX (1) = ' , F9 . 6 , / , 6X, ' EY (1) = ' , F9.6, / , 

6X, 'EX (2) = ’ ,F9.6,/,6X, 'EY (2) = ' ,F9.6,/, 

6X, 'CALL CURVE(EX,EY, 2,0)') 

8000 FORMAT (6X, ' XO ( ' , 12 , ' ) = 1 , F9 . 6 , / , 6X, ’ YO ( ' , 12, ' )= ' , F15 . 6) 

8001 FORMAT (6X, 'XI (' ,12, ' ) = ' , F9 . 6 , / ,6X, ' Y1 ( ' , 12, ’)=',F15.6) 

9000 FORMAT (6X, 'XOC ,12, ')=' ,F7.3,/,6X, ' YO ( ’ ,12, ')=' , F16 . 4) 

9001 FORMAT (6X, ’ XI ( 1 , 12 , ' ) = ' , F7 . 3 , /, 6X, ' Yl ( ’ , 12, ' ) = ’ , F16 . 4) 

END 

A 

* for THE DETERMINATION OF MEAN CONTACT POINT 

5 1 : 

SUBROUTINE PCN1(X,F,NE) 

IMPLICIT REAL*8(A-H,K,M-Z) 

CHARACTER *8 HG 
INTEGER NE 

REAL *8 X (NE) , F (NE) , PAR (6) 

COMMON /P 1 /PAR 
COMMON/AO/HG 

COMMON/ Al/pll,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3 
COMMON /a2 /SG, CSRT2 , QG, SNPSIG, CSPSIG 
COMMON/A3/TND1 , TND2 , RITAG 
THETAG=X(1) 

CSTHEG=DCOS (THETAG) 

SNTHEG=DS IN (THETAG) 

IF (HG .EQ. ' L ' ) THEN 

UG=PAR (1) -SG* (PAR (2) *SNTHEG+PAR (3) ) / (CSRT2*DSIN (THETAG-QG) ) 
ELSE 

UG=PAR (1) -SG* (PAR (2) *SNTHEG+PAR (3) ) / (CSF T 2*DSIN (THETAG+QG) ) 
END IF 

Bcx=PAR (4)-UG*CSPSIG 
Bcy=UG*SNPSIG*SNTHEG 
Bcz=UG*SNPSIG*CSTHEG 
CALL TRCOOR (Bpx , Bpy , Bpz, 

. pll,pl2,pl3,p21,p22, p23 , p31 , p32 , p33, pi ,p2,p3, 

. Bcx,Bcy,Bcz) 

XM=Bpz* (TND1-TND2) /2 . DOO 

F(l)=Bpx-XM 

END 

* 

* FOR THE DETERMINATION OF MEAN CONTACT POINT 

* 

SUBROUTINE PCN2(X,F,NE) 

IMPLICIT REAL*8(A-H,K,M-Z) 

CHARACTER *8 HG 
INTEGER NE 

REAL*8 X (NE) , F (NE) , PAR (6) 

COMMON/P 1 /PAR 
COMMON/AO/HG 

COMMON/Al/pll , p 12,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3 
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C - 



C0MM0N/A2/SG , CSRT2 , QG , SNPSIG , CSPSIG 
COMMON / A3 / TND 1 , TND2 , R I TAG 
COMMON /A4/CSD2 , SND2 , CSPIT2 , SNPIT2 
COMMON /A5/CSQG, SNQG, THETAG 
PHIG=X (1) 

CSPHIG=DCOS (PHIG) 

SNPHIG=DSIN (PHIG) 

IF (HG .EQ. ' L ' ) THEN 

UG=PAR (l)-SG* (PAR (2) +PAR (3) *DSIN (QG-PHIG) ) / 

(CSRT2’ V DS IN (THETAG-QG+PH I G) ) 

ELSE 

UG=PAR (1) -SG* (PAR (2) +PAR (3) *DSIN (QG-PHIG) ) / 

(CSRT2*DSIN (THETAG+QG-PHIG) ) 

END IF 

Bcx=PAR(4)-UG*CSPSIG 
Bcy=UG*PAR (5) 

Bcz=UG*PAR (6) 

* 

* Mmc=Mms*Ms c 

ix 

IF (HG .EQ. ' L 1 ) THEN 

CALL COMBI (mll,ml2,ml3,m21 , m22 ,m23 ,m31 ,m32 ,m33 ,ml ,m2,m3 , 

1 . DOO , 0. DOO , 0. DOO , 0 .DOO , CSPHIG, SNPHIG, O.DOO, -SNPHIG, CSPHIG, 

. O.DOO, O.DOO, O.DOO, 

1 . DOO , 0 . DOO , 0 . DOO , 0 . DOO , CSQG , -SNQG , 0 . DOO , SNQG , CSQG , 

. O.DOO, -SG*SNQG,SG*CSQG) 

Vc 

* Mpc-Mpm*Minc 

CALL COMBI (pll,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3, 

CSD2 , 0 . DOO , -SND2 , 0 . DOO , 1 . DOO , 0 . DOO , SND2 , 0 . DOO , CSD2 , 

. O.DOO, O.DOO, O.DOO, 

. mll,ml2,ml3,m21 ,m22 ,m23 , m31 ,m32 ,m33 , ml , m2 , m3) 

ELSE 

* 

* Mmc=Mms*Ms c 

Vc 

CALL COMB I (m 1 1 , m 1 2 , m 1 3 , m2 1 , m22 , m23 , m3 1 , m32 , m33 , m 1 , m2 , m3 , 

1. DOO, O.DOO, O.DOO, O.DOO, CSPHIG, -SNPHIG, O.DOO, SNPHIG, CSPHIG, 

. O.DOO, O.DOO, O.DOO, 

1. DOO, O.DOO, O.DOO, O.DOO, CSQG, SNQG, O.DOO, -SNQG, CSQG, 

. 0 . DOO , SG*SNQG , SG*CSQG) 

* 

* Mpc = MpTn ,f Mmc 

Vc 

CALL COMBI (pll,pl2,pl3,p21 ,p22,p23,p31 ,p32,p33,pl,p2,p3, 
CSD2 , 0 . DOO , -SND2 , 0 . DOO , 1 . DOO , 0 . DOO , SND2 , 0 . DOO , CSD2 , 

. O.DOO, O.DOO, O.DOO, 

. mll,ml2,ml3, m2 I,m22,m23,m31,m32, m3 3 , ml , m2 , m3) 

END IF 

CALL TRCOOR (Bpx, Bpy , Bpz, 

. pll,pl2,pl3,p21 ,p22,p23,p31 ,p32,p33,pl,p2,p3, 

. Bcx,Bcy,Bcz) 
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XM=Bpz* (TND1 -TND2) /2 . DOO 

F(l)=Bpx-XM 

RETURN 

END 

* 

* F0R THE DETERMINATION of coordinates and normals of contact points 


SUBROUTINE TCN (TX , TF , NE) 

IMPLICIT REAL *8 (A~H , K , M~2) 

REAL*8 LANDAP , LM 
CHARACTER *8 HG 
INTEGER NE 

DIMENSION TX(NE) ,TF(NE) ,TPAR(19) 

COMMON/T1/TPAR 

COMMON/AO/HG 

COMMON/ A2/ SG , CSRT2 , QG , SNPS IG , CSPS IG 
COMMON/ A4/ CSD2 , SND2 , CSPIT2 , SNPIT2 

COMMON/B 1 /CSPH1 1 , SNPH1 1 , SP , EM , LM , CSRT1 , CSD1 , SND1 , CSLANP , SNLANP 
COMMON/B2/CSPIT1 , SNPITl , MP1 ,MG2,QP 
COMMON/B3/B2f x , B2f y , B2f z 

COMMON/B4/CSPH2 , SNPH2 , CSPH2 1 , SNPH21 
COMMON/B5/XCR ,ZCR 
COMMON/C 1 /UG, CSTAUG, SNTAUG 
COMMON/C2/N2fx,N2fy,N2£z 

COMMON/D1/CSTAUP , SNTAUP 
COMMON/F1/PHIGO 
COMMON/G1/DA, DV 
J=IDINT (TPAR (16) ) 

PHIP=TX (1) 

THETAP=TX (2) 

LANDAP=TX (3) 

PHI21 = TX (4) 

PHIG=TX (5) 

THETAG=TX (6) 

CSPHIP=DCOS (PHIP) 

SNPHIP=DSIN (PHIP) 

CSTHEP=DCOS (THETAP) 

SNTHEP=DSIN (THETAP) 

CSLANP=DCOS (LANDAP) 

SNLANP=DSIN (LANDAP) 

CSPH21=DCOS (PHI21) 

SNPH21=DSIN (PHI21) 

CSPHIG=DCOS (PHIG) 

SNPHIG=DSIN (PHIG) 

CSTHEG=DCOS (THETAG) 

SNTHEG=DSIN (THETAG) 

PHI2= (PHIG-PHIGO) /MG2 
PHI 1=PHIP/MP1 
CSPH2=DCOS (PHI2) 

SNPH2=DSIN(PHI2) 

CSPHl=DCOS (PHI1) 

SNPH1=DSIN (PHI1) 
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IF (HG .EQ. 1 L ' ) THEN 
TAUP=THETAP+QP-PHIP 
ELSE 

TAUP=THETAP-QP+PHIP 
END IF 

CSTAUP=DCOS (TAUP) 

SNTAUP=DSIN (TAUP) 

IF (HG .EQ. ' L ' ) THEN 
TAUG=THETAG-QG+PHIG 
ELSE 

TAUG= THETAG+ QG-PH I G 
END IF 

CSTAUG=DCOS (TAUG) 

SNTAUG=DSIN (TAUG) 

CSQPHP=DCOS (QP-PHIP) 

SNQPHP=DSIN (QP-PHIP) 

CSQPHG=DCOS (QG-PHIG) 

SNQPHG=DSIN (QG-PHIG) 

Vc 

* LEFT-HAND GEAR 

* 

IF (HG .EQ. ' L 1 ) THEN 

UG=TPAR (1) -SG* (TPAR (2) *SNTHEG~SNQPHG*TPAR (3) ) / (CSRT2*SNTAUG) 
B2py=UG*SNPSIG*SNTAUG-SG*SNQPHG 

ELSE 

UG=TPAR (1) -SG* (TPAR (2) *SNTHEG+SNQPHG*TPAR (3) ) / (CSRT2*SNTAUG) 
B2py=UG*SNPSIG*SNTAUG+SG*SNQPHG 

END IF 

B2px=CSD2* (TPAR (4) -UG*CSPSIG) -SND2* (UG*SNPSIG*CSTAUG+SG’ r CSQPHG) 
B2pz=SND2* (TPAR (4) -UG*CSPSIG) +CSD2* (UG*SNPSIG*CSTAUG+SG*CSQPHG) 
N2px=TPAR (5) -TPAR (6) *CSTAUG 
N2py=CSPS IG*SNTAUG 
N2pz=TPAR (7) +TPAR (8) *CSTAUG 

* 

* [Mwp] = [Mwa] [Map] 

* 

IF (HG .EQ. ' L ' ) THEN 

CALL COMBI (wpll,wpl2,wpl3,wp21,wp22,wp23,wp31,wp32,wp33, 

« wp 1 , wp2 , wp3 , 

*. CSPH2 , SNPH2 ’ 0 . D00 , -SNPH2 , CSPH2 , 0 . D00 , 0 . DOO , 0 . DOO , 1 . DOO , 

. 0. DOO, 0. DOO, 0. DOO, 

. CSPIT2,0.D00,SNPIT2,0.D00, 1 .DOO,O.DOO,-SNPIT2,O.DOO,CSPIT2, 

. 0. DOO, 0. DOO, 0. DOO) 

ELSE 

CALL COMB I (wp 1 1 , wp 1 2 , wp 1 3 , wp2 1 , wp22 , wp23 , wp3 1 , wp32 , wp33 , 
wpl ,wp2,wp3, 

. CSPH2 , -SNPH2 , 0 . DOO , SNPH2 , CSPH2 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

. 0. DOO, 0. DOO, 0. DOO, 

. CSPIT2,0.D00,SNPIT2,0.D00, 1 .D00,0.D00,-SNPIT2,0.D00,CSPIT2, 

. 0. DOO, 0. DOO, 0. DOO) 


187 


END IF 

CALL TRCOOR (B2wx , B2wy , B2wz, 

. wp 1 1 , wp 12 , wp 13 , wp2 1 , wp22 , wp23 , wp31 , wp32 , wp33 , wp 1 , wp2 , wp3 , 

. B2px , B2py , B2pz) 

CALL TRCOOR (N2wx , N2wy , N2wz , 

. wpl 1 , w P 12 , w P 13 , w P 21 , w P 22 , w P 23 , w P 31 , w P 32, w P 33 , 0 . D00 , 0 . D00 , 0 . D00, 
. N2px , N2py , N2pz) 

[Mf w] = [Mf a] [Maw] 

IF (HG .EQ. ' L ' ) THEN 

CALL COMB I (fwl 1 , f wl2 , f wl3 , f w2 1 , f w22 , f w23 , fw3 1 , f w32 , f w33 , 

. fwl,fw2,fw3, 

. CSPIT2 , 0. D00 , -SNPIT2 , 0. D00 , 1 . DOO, O.DOO, SNPIT2, 0. DOO, CSPIT2 , 

. 0 . D00 , 0 . D00 , 0 . DOO , 

. CSPH21 , -SNPH21 , 0. DOO , SNPH21 , CSPH21 ,0.D00, O.DOO , O.DOO 1 .DOO 
. O.DOO, O.DOO, O.DOO) 

ELSE 

CALL COMBI (f wl 1 , f wl2 , fwl3 , f w21 , f w22 , f w23 , f w31 , f w32 , fw33 , 

• fwl,fw2,fw3, 

. CSPIT2, 0. DOO, -SNPIT2, O.DOO, l.DOO,O.DOO,SNPIT2,O.DOO,CSPIT2, 

. O.DOO, O.DOO, O.DOO, 

. CSPH21 ,SNPH21 , O.DOO, -SNPH21 ,CSPH21, O.DOO, O.DOO, O.DOO, 1 .DOO, 

. O.DOO, O.DOO, O.DOO) 

END IF 

CALL TRCOOR (B2fx,B2fy,B2fz, 

. f wl 1 , f wl2 , f wl 3 , f w2 1 , f w22 , f w23 , f w3 1 , f w32 , f w33 , f wl , f w2 , f w3 , 

. B2wx , B2wy , B2wz) 

CALL TRCOOR (N2fx,N2fy,N2fz, 

. fwl I,fwl2,fwl3,fw21,fw22,fw23,fw31,fw32,fw33, O.DOO, O.DOO, O.DOO, 

. N2wx,N2wy,N2wz) 

PINION 

IF (HG .EQ. ' L ' ) THEN 

Blpy= (ZCR+TPAR (17) *SNLANP) *SNTAUP+SP*SNQPHP-EM 
ELSE 

Blpy= (ZCR+TPAR (17) *SNLANP) ,V SNTAUP-SP*SNQPHP+EM 
END IF 

Blpx= (XCR+TPAR (17) *CSLANP) *CSD1- ( (ZCR+TPAR (17) *SNLANP) *CSTAUP+ 
SP*CSQPHP+LM)*SND1 

Blpz= (XCR+TPAR (17) *CSLANP) *SND1+ ( (ZCR+TPAR (17) *SNLANP) *CSTAUP+ 
SP*CSQPHP+LM) *CSD1 
IF ( J .EQ. 2) THEN 

Nlpx=CSLANP*CSDl-SNLANP*SNDl ,v CSTAUP 
N 1 py=SNLANP *SNTAUP 
Nlpz=CSLANP*SNDl+SNLANP*CSDl*CSTAUP 
ELSE 

Nlpx=-CSLANP*CSD1+SNLANP*SND1*CSTAUP 
Nlpy=-SNLANP*SNTAUP 
Nlpz=-CSLANP ,V SND1-SNLANP*CSD1*CSTAUP 
END IF 


* [Mwp] = [Mwa] [Map] 

Vc 

IF (HG .EQ. ' L ' ) THEN 

CALL COMB I (wp 1 1 , wp 1 2 , wp 1 3 , wp2 1 , wp22 , wp23 , wp3 1 , wp32 , wp33 , 

, wp 1 j wp 2 y wp 3 , 

. CSPHl , -SNPHl , 0. D00, SNPHl , CSPHl , 0.D00 , O.DOO , O.DOO, 1 .D00 , 

. O.DOO, O.DOO, O.DOO, 

. CSPITl , O.DOO, SNPITl, O.DOO, 1 .D00, O.DOO, -SNPIT1 , O.DOO, CSPITl , 

. O.DOO, O.DOO, O.DOO) 

ELSE 

CALL COMBI (wpll , wpl2 , wpl3 , wp21 , wp22 , wp23, wp31 , wp32 , wp33 , 

„ wpl,wp2,wp3, 

. CSPHl , SNPHl , O.DOO, -SNPHl , CSPHl, O.DOO, O.DOO, O.DOO, 1 .DOO, 

. O.DOO, O.DOO, O.DOO, 

. CSPITl , O.DOO, SNPITl .O.DOO, 1 .DOO, O.DOO, -SNPITl , 0 .DOO, CSPITl , 

. O.DOO, O.DOO, O.DOO) 

END IF 

CALL TRCOOR (Blwx, Blwy, Blwz, 

. wpl 1 , w P 12 , wpl3 , wp21 , w P 22 , w P 23 , w P 31 , wp32, wp33 , wpl , w P 2, wp3 , 

. Blpx,Blpy,Blpz) 

CALL TRCOOR (Nlwx , Nlwy , Nlwz, 

. wp 1 1 , wp 12 , wp 1 3 , wp2 1 , wp22 , wp23 , wp3 1 , wp32 , wp33 , 0 . DOO , 0 . DOO , 0 . DOO , 

. Nlpx.Nlpy.Nlpz) 

* 

* [Mpw] = [Mpa] [Maw] 

* 

IF (HG .EQ. ' L ' ) THEN 

CALL COMB I (pwl 1 , pwl 2 , pwl 3 , pw2 1 , pw22 , pw23 , pw3 1 , pw32 , pw33 , 

• pwl , pw2 y pw3 y 

. CSPITl , O.DOO, -SNPITl ,0.D00, 1 . DOO , 0 . DOO, SNPITl , O.DOO , CSPITl , 

. O.DOO, O.DOO, O.DOO, 

. CSPHl 1 , SNPH 11,0. DOO , -SNPHl 1 , CSPHl 1 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 

. O.DOO, O.DOO, O.DOO) 

ELSE 

CALL COMBI (pwl 1 , pwl2 , pwl3 , pw21 , pw22, pw23 , pw31 , pw32, pw33 , 

• pwl y pw2 y pw3 y 

. CSPITl , O.DOO, -SNPITl , O.DOO, 1 . DOO , 0. DOO, SNPITl , O.DOO, CSPITl , 

. O.DOO, O.DOO, O.DOO, 

. CSPHl 1, -SNPHl 1,0. DOO, SNPHl 1, CSPHl 1,0. DOO, 0. DOO, 0. DOO, 1. DOO, 

. O.DOO, O.DOO, O.DOO) 

END IF 

CALL TRCOOR (Blpx , B lpy , Blpz, 

. pwl 1 , pwl2 , pwl3 , pw2 1 , pw22 , pw23 , pw3 1 , pw32 , pw33 , pwl , pw2 , pw3 , 

. Blwx, Blwy, Blwz) 

CALL TRCOOR (Nlpx.Nlpy.Nlpz, 

. pwl 1 , pwl 2 , pwl 3 , pw21 , pw22 , pw23 , pw3 1 , pw32 , pw33 , 0 . DOO, 0 . DOO , 0 . DOO , 
. Nlwx, Nlwy , Nlwz) 

Blfx=-Blpx+DA*SNPIT1 

Blfy=-Blpy+DV 

Blfz=Blpz+DA*CSPITl 

Nlfx=-Nlpx 

Nlfy=-Nlpy 

Nlfz-Nlpz 
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IF(HG .EQ. ' L 1 ) THEN 

TF (1) = (TPAR (9) *SNTAUP+TPAR (10) *SNQPHP~TPAR (11)) *CSLANP~ 

(TPAR (12) *SNTAUP~TPAR (13) *SNTHEP+TPAR (14) *CSTAUP+TPAR (15) * 
SNTAUP) *SNLANP 

ELSE 

TF (1) = (TPAR (9) *SNTAUP-TPAR (10) *SNQPHP+TPAR (1 1) ) *CSLANP- 

(TPAR (12) *SNTAUP-TPAR (13) *SNTHEP~TPAR (14) *CSTAUP+TPAR (15) * 
SNTAUP) *SNLANP 

END IF 

TF (2) =B2f x-Blf x 
TF(3)=B2fy-Blfy 
TF(4)=B2fz-Blfz 
TF(5)=N2fx-Nlfx 
TF(6)=N2fy-Nlfy 
END 

* FOR THE DETERMINATION OF GEAR PRINCIPAL CURVATURES AND DIRECTIONS 

* 

SUBROUTINE PRING2 (KSG, GG, EGIfx , EGIfy , EGIf z, EGI If x , EGIIfy, EGIIf z) 
IMPLICIT REAL ,V 8(A-H,K,M-Z) 

COMMON/ A2/SG, CSRT2 , QG , SNPSIG, CSPSIG 
COMMON/A3/TND1 , TND2 .RITAG 
COMMON/A4/CSD2 , SND2 , CSPIT2 , SNPIT2 
COMMON/B2/CSPIT1 , SNPIT1 ,MP1 ,MG2,QP 
COMMON/B3/B2fx,B2fy,B2fz 
COMMON/B4/ CSPH2 , SNPH2 , CSPH21 , SNPH21 
COMMON/C 1 /UG , CSTAUG , SNTAUG 
COMMON/C2/N2fx,N2fy,N2fz 
COMMON/ FI /PH I GO 
KCI=-CSPSIG/ (UG*SNPSIG) 

KCII=0.D00 

ECIfx=SND2*SNTAUG 

ECIfy=CSTAUG 

ECIfz=-CSD2*SNTAUG 

ECIIfx=-CSD2*CSPSIG-SND2*SNPSIG*CSTAUG 

ECIIfy=SNPSIG*SNTAUG 

ECI If z =- SND2 ,c CSPSIG+CSD2 ,v SNPSIG*CSTAUG 

WGfx=-SNPIT2 

WGfy=0.D00 

WGfz=CSPIT2 

WCfx=-MG2*CSD2 

WCfy=0.D00 

WCfz=-MG2*SND2 

WGCf x=WGf x-WCf x 

WGCfy=WGfy-WCfy 

WGCfz=WGfz-WCfz 

CALL CROSS (VTGf x , VTGf y , VTGf z, WGf x , WGfy , WGf z, B2f x , B2f y , B2fz) 

CALL CROSS (VTCfx , VTCfy , VTCfz, WCf x , WCfy , WCfz,B2fx , B2fy , B2fz) 

VTGCf x=VTGf x-VTCfx 
VTGCfy=VTGfy- VTCfy 
VTGCf z=VTGfz-VTCfz 

CALL DOT (VCI , EC I f x , EC I f y , ECI f z , VTGCf x , VTGCf y , VTGCf z) 

CALL DOT (VCI I, ECI I fx.ECIIfy, ECI If z, VTGCf x, VTGCf y, VTGCf z) 


190 


* CALCULATE Al3,A23,A33 

* 

CALL DET (DETI , WGCfx, WGCf y, WGCf z,N2fx , N2fy , N2fz, ECIfx, ECIfy , ECIfz) 
A13=-KCI*VCI-DETI 

CALL DET (DETII, WGCfx, WGCf y, WGCf z,N2fx,N2fy,N2fz, 
ECIIfx,ECIIfy,ECIIfz) 

A23=-KC I I *VC I I-DET I I 

CALL DET (DET3, N2f x,N2fy,N2fz, WGCfx, WGCf y, WGCf z, 

VTGCf x , VTGCfy , VTGCfz) 

CALL CROSS (Cx , Cy , Cz , WGf x , 0 . D00 , WGf z , VTCf x , VTCf y, VTCf z) 

CALL CROSS (Dx , Dy , Dz, WCf x , 0 . DOO , WCfz, VTGf x , VTGfy, VTGf z) 

CALL DOT (DET45 , N2 f x , N2f y , N2 f z , Cx-Dx , Cy-Dy , Cz~Dz) 
A33=KCI*VCI*VCI+KCII*VCI I ,r VCI I-DET3-DET45 

* 

* CALCULATE SIGMA 

jV 

P=A23*A23~A13*A13+ (KCI-KCI I) *A33 
SIGMA2=DATAN (2 . D00*A13*A23/P) 

SIGMA=0.5D00*SIGMA2 

* 

* CALCULATE KG I AND KGII 

'k 

GG=P/ (A33*DCOS (SIGMA2) ) 

KSG=KC I +KC 1 1 - ( A 1 3 * A 1 3+ A23 * A2 3 ) / A3 3 
KGI= (KSG+GG) /2, DOO 
KGI 1= (KSG-GG) / 2 . DOO 

* 

* CALCULATE EGI AND EGII 

* 

CALL ROTATE (EGI f x, EGI f y, EGI fz, ECIfx, ECIfy, ECIfz, -SIGMA, N2fx,N2fy, 
. N2fz) 

CALL ROTATE(EGIIfx,EGIIfy,EGIIfz,EGIfx,EGIfy,EGIfz,RITAG, 

. N2fx,N2fy,N2fz) 

END 

* 

* FOR THE DETERMINATION OF PINION PRINCIPAL CURVATURES AND DIRECTIONS 

* 

SUBROUTINE PRINP1 (KSP , GP , EPIf x , EPIf y , EPIfz, EPI If x , EPII f y , EPI If z) 
IMPLICIT REAL*8(A-H,K,M-Z) 

REAL*8 LM, TPAR (19) 

COMMON / T 1 / TP AR 

COMMON / A3 / TND 1 , TND2 , R I TAG 

COMMON/B 1 /CSPH1 1 , SNPHl 1 , SP , EM , LM , CSRT1 , CSD1 , SND1 , CSLANP , SNLANP 

C0MM0N/B2/CSPIT1 , SNPIT1 , MP 1 ,MG2 , QP 

COMMON/B 3 /B2fx , B2fy , B2f z 

COMMON/B5/XCR , ZCR 

COMMON/C2/N2fx,N2fy,N2fz 

COMMON/Dl / CSTAUP , SNTAUP 

J=IDINT (TPAR (16) ) 

R=TPAR (17) 

IF (J .EQ. 1)THEN 
KCI=SNLANP/ (ZCR+R*SNLANP) 
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KCII=1 .DOO/R 
ELSE 

KCI=-SNLANP/ (ZCR+R*SNLANP) 

KCII=-1. DOO/R 
END IF 

ECI f x=SNDl*SNTAUP 

ECIfy=CSTAUP 

ECIfz-CSDl*SNTAUP 

EC I If x=CSDl*SNLANP+SNDl*CSLANP*CSTAUP 
EC I I f y=-CSLANP*SNTAUP 

EC 1 I f z=-SNDl *SNLANP+CSDl *CSLANP*CSTAUP 

IF ( J .EQ. 2) THEN 

ECIIf x=~ECI If x 

ECIIfy=-ECIIfy 

ECIIf z= _ ECI If z 

END IF 

WPfx=-SNPITl 

WPfy=O.DOO 

WPfz=-CSPITl 

WCfx=-MPl*CSDl 

WCfy=O.DOO 

WCfz=MPl 5r SNDl 

WPCf x=WPf x-WCf x 

WPCfy=WPfy-WCfy 

WPCfz=WPfz-WCfz 

CALL CROSS (VTPf x, VTPfy, VTPf z, WPf x , WPf y , WPf z, B2f x , B2f y , B2f z) 

CALL CROSS (VTC1 f x , VTClfy , VTClf z, WCf x , WCfy, WCf z, B2f x, B2fy , B2f z) 

CALL CROSS (VTC2f x , VTC2fy , VTC2f z, LM*SNDl , EM, LM*CSD1 , WCf x , WCf y, WCf z) 

VTCf x=VTCl f x+VTC2f x 

VTCfy=VTClfy+VTC2fy 

VTCfz=VTClfz+VTC2fz 

VTPCf x=VTPf x-VTCf x 

VTPCf y= VTPfy- VTCfy 

VTPCf z=VTPfz-VTCfz 

CALL DOT (VC I ,ECIfx,ECIfy,ECIfz, VTPCf x, VTPCf y, VTPCf z) 

CALL DOT (VCI I , ECI I f x , ECI If y , ECI I fz, VTPCf x , VTPCf y , VTPCf z) 

* 

* CALCULATE Al3,A23,A33 

* 

CALL DET (DETI , WPCf x , WPCfy , WPCf z,N2f x,N2fy ,N2fz, ECIf x, ECIfy , ECIf z) 
A13=-KCI*VCI-DETI 

CALL DET (DETI I, WPCf x, WPCfy, WPCf z,N2fx,N2fy,N2fz, 
ECIIfx.ECIIfy.ECIIfz) 

A23=-KCII*VCII-DETII 

CALL DET (DET3 , N2f x , N2f y , N2f z , WPCf x , WPCfy , WPCf z, 

VTPCf x , VTPCf y, VTPCf z) 

CALL CROSS (Cx , Cy , Cz, WPf x , 0 . DOO , WPf z, VTCf x , VTCfy , VTCf z) 

CALL CROSS (Dx , Dy , Dz , WCf x , 0 . DOO , WCf z, VTPf x .VTPfy , VTPf z) 

CALL DOT (DET45 , N2f x , N2f y , N2f z , Cx~Dx , Cy-Dy , Cz~Dz) 
A33=KCI*VCI*VCI+KCII*VCII*VCII-DET3-DET45 
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* CALCULATE SIGMA 

* 


P=A23*A23-A13*A13+ (KCI-KCII) *A33 
SIGMA2=DATAN (2 . D00*Al3*A23/P) 

S I GMA= 0 . 5 DOO * S I GMA2 

'V 

* CALCULATE KPI AND KPII 

* 

GP=P/ (A33*DCOS (SIGMA2) ) 

KSP=KCI+KCII-(A13*A13+A23*A23)/A33 
KPI= (KSP+GP) /2 . DOO 
KPII= (KSP-GP) /2.D00 

* 

* CALCULATE EPI AND EPII 

* 

CALL ROTATE (EPIf x, EPIfy, EPIf z,ECIfx,ECIfy,ECIfz, “SIGMA, N2fx,N2fy, 
. N2fz) 

CALL ROTATE (EPII fx, EPI Ify,EPIIfz, EPIf x , EPIfy , EPIf z, RITAG , 

. N2fx,N2fy,N2fz) 

END 

it 

* FOR THE DETERMINATION OF THE ANGLE BETWEEN GEAR PRINCIPAL DIRECTIONS 

* AND PINION PRINCIPAL DIRECTIONS 

SUBROUTINE SIGAN2 (EGIf x , EGI f y , EGIf z, EGIIf x , EGI Ify , EGIIf z, EPI f x , 

. EPIfy, EPIf z,CS2SIG,SN2SIG,SIGMPG) 

IMPLICIT REAL ,V 8(A-H,K,M-Z) 

CALL DOT (CSSIG , EPI f x , EP I f y , EP I f z , EGI f x , EGI f y , EGIf z) 

CALL DOT (SNSIG,EPIfx,EPIfy,EPIfz, - EGIIfx, -EGI Ify , -EGIIf z) 

SIGM2=4 .D00*DATAN (SNSIG/ (1 .DOO+CSSIG) ) 

S IGMPG= . 5D00*S I GM2 
CS2SIG=DCOS (SIGM2) 

SN2SIG=DSIN (SIGM2) 

END 

it 

* FOR THE DETERMINATION OF CONTACT ELLIPS 

it 

SUBROUTINE ELLIPS (KSG , GG, KSP , GP , CS2SIG, SN2SIG, DEF .ALFAP , 

AXISL, AXISS, EPIf x, EPIf y, EPIf z) 

IMPLICIT REAL*8 (A-H,K,M-Z) 

COMMON/ A3 /TND1 , TND2 , RITAG 
COMMON/C2/N2fx,N2fy,N2fz 
COMMON/E 1 /XB f , YBf , ZBf 

D=DSQRT (GP*GP-2 . D00*GP*GG*CS2SIG+GG*GG) 

CS2AFP= (GP-GG*CS2SIG) /D 
SN2AFP=GG*SN2SIG/D 
ALFAP=DATAN(SN2AFP/ (1 . D00+CS2AFP) ) 

A= . 25D00*DABS (KSP-KSG-D) 

B= . 25D00*DABS (KSP-KSG+D) 

IF (KSG .LT. KSP) THEN 
AXISL=DSQRT (DEF/A) 

AXISS=DSQRT (DEF/B) 

CALL ROTATE (XBf , YBf , ZBf ,EPIfx, EPIfy, EPI fz , RITAG- ALFAP , N2fx , 
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. N2fy,N2fz) 

ELSE 

AXISL=DSQRT (DEF/B) 

AXISS=DSQRT (DEF/A) 

CALL ROTATE (XBf , YBf , ZBf , EPIf x , EPIfy , EPIf Z, -ALFAP , N2f x , N2f y , 

. N2fz) 

END IF 

XBf =AXISL ,v XBf 
YBf =AXISL*YBf 
ZBf =AXISL*ZBf 
END 

* 

* COORDINATE TRANSFORMATION FOR F TO P 

* 

SUBROUTINE PF (B2px , B2py , B2pz, B2f x , B2fy , B2f z) 

IMPLICIT REAL*8(A-H,K,M-Z) 

COMMON /A4/CSD2 , SND2 , CSPIT2 , SNPIT2 
COMMON/B4/CSPH2 , SNPH2.CSPH21 , SNPH21 

* 

* [Mtf] = [Mta] [Maf] 

Vc 

CALL COMBI (tll,tl2,tl3,t21,t22,t23,t31,t32,t33,tl,t2,t3, 

. CSPH2 1 , SNPH2 1 , 0 . DOO , -SNPH2 1 , CSPH2 1 , 0 . DOO , 0 . DOO , 0 . DOO , 1 . DOO , 
. 0. DOO, 0. DOO, 0. DOO, 

. CSPIT2 , 0 . DOO , SNPIT2 , 0 . DOO, 1 .DOO,O.DOO,-SNPIT2,O.DOO,CSPIT2, 
. 0. DOO, 0. DOO, 0. DOO) 

CALL TRCOOR (B2wx , B2wy , B2wz, 

. tll,tl2,tl3,t21,t22,t23,t31,t32,t33,tl,t2,t3, 

. B2f x , B2fy , B2f z) 

* 

* [Mpt] = [Mpa] [Mpt] 

* 

CALL COMBI (pll,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3, 

. CSP IT2 , 0 . DOO , -SNP I T2 , 0 . DOO , 1 . DOO , 0 . DOO , SNP I T2 , 0 . DOO , CSP IT2 , 
. 0. DOO, 0. DOO, 0. DOO, 

. CSPH2 , -SNPH2 , 0 . DOO , SNPH2 , CSPH2, O.DOO.O.DOO, 0. DOO , 1 . DOO , 

. 0. DOO, 0. DOO, 0. DOO) 

CALL TRCOOR (B2px,B2py,B2pz, 

. p ll,pl2,pl3,p21,p22,p23,p31,p32,p33,pl,p2,p3, 

. B2wx , B2wy , B2wz) 

END 

* 

* USING EULER FORMULA TO DETERMINATION SURFACE INTERFERENCE 

* 

SUBROUTINE EULER (KSG , GG , KSP , GP , CS2SIG, SN2SIG , IEU) 

IMPLICIT REAL*8(A-H,K,M-Z) 

A=KSG-KSP 

B=DSQRT ( (GG-GP*CS2SIG) **2+ (GP*SN2SIG) **2) 

KR1= (A+B) /2.D00 

KR2=(A"B)/2.D00 

IF (KRl*KR2 .LT. 0. DOO) THEN 

IEU=1 

ELSE 
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IEU=0 
END IF 
END 


* DETERMINANT 

* 

SUBROUTINE DET (S , A, B , C , D, E , F , G, H, P) 

IMPLICIT REAL*8(A-H,K,M-Z) 

S=A*E*P+D*H*C+G*B*F-A*H*F-D*B*P-G*E*C 

RETURN 

END 

* 

* COORDINATE TRANSFORMATION 

* 

SUBROUTINE TRCOOR (XN , YN , ZN , R1 1 , R12 , R13 , R21 , R22 ,R23 , R31 , R32 ,R33 , 
T1,T2,T3,XP,YP,ZP) 

IMPLICIT REAL*8(A-H,0-Z) 

XN=R11*XP+R12*YP+R13*ZP+T1 

YN=R2 1 *XP+R22*YP+R23*ZP+T2 

ZN=R31*XP+R32*YP+R33*ZP+T3 

RETURN 

END 

* 

* MULTIPLICATION OF TWO TRANSFORMATION MATRICES 

* 

SUBROUTINE COMBI (Cl 1 , C12 , C13 , C21 , C22 , C23 ,C31 , C32 , C33 , Cl ,C2 , C3 , 
All .AiZ.An.AZl.AZZ.AZS.ASl.ASZ.ASS.Al.AZ.AS, 
B11,B12,B13,B21, B22 , B23 , B31 , B32,B33 , B1 , B2, B3) 
IMPLICIT REAL*8(A-H,M,N,0-Z) 

C11=B31*A13+B21*A12+B11*A11 

C12=B32*A13+B22*A12+B12*A11 

C13=B33*Al3+B23*Al2+Bl3*All 

C21=B31*A23+B21*A22+B11*A21 

C22=B32*A23+B22*A22+B 1 2*A2 1 

C23=B33’ v A23+B23*A22+B13*A21 

C3 1 =B3 1 *A33+B2 l *A32+B 1 l *A3 1 

C32=B32*A33+B22*A32+B12*A31 

C33=B33*A33+B23*A32+B13*A31 

C1=B3*A13+B2*A12+B1*A11+A1 

C2=B3*A23+B2*A22+B1*A21+A2 

C3=B3*A33+B2*A32+B 1 *A3 1+A3 

RETURN 

END 

* 

* DOT OF TWO VECTORS 

* 

SUBROUTINE DOT (V , Xl , Yl , Z 1 , X2 , Y2 , Z2) 

IMPLICIT REAL*8(A-H,0-Z) 

V=X1*X2+Y1*Y2+Z1*Z2 

RETURN 

END 

* 

* CROSS OF TWO VECTORS 
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SUBROUT INE CROSS (X,Y,Z,A,B,C,D,E,F) 

IMPLICIT REAL*8(A-H,0-Z) 

X=B*F~C*E 

Y-C*D-A*F 

Z=A*E~B*D 

RETURN 

END 

ROTATION A VECTOR ABOUT ANOTHER VECTOR 

SUBROUTINE ROTATE (XN , YN , ZN , XP , YP , ZP , THETA, UX , UY , UZ) 

IMPLICIT REAL*8(A-H,0-Z) 

CT=DCOS (THETA) 

ST=DSIN (THETA) 

VT=1.DOO-CT 

Rl 1=UX*UX*VT+CT 

R12=UX*UY*VT-UZ*ST 

R13'UX*UZ*VT+UY*ST 

R21=UX*UY*VT+UZ*ST 

R22=UY*UY*VT+CT 

R23-UY*UZ*VT~UX*ST 

R31=UX*UZ*VT-UY*ST 

R32«UY*UZ*VT+UX*ST 

R33=UZ*UZ*VT+CT 

CALL TRCOOR (XN, YN , ZN , Rl 1 ,R12,R13,R21 ,R22,R23 ,R31 ,R32 , R33 , 

0 . D00 , 0 . DOO , 0 . D00 , 

XP , YP , ZP) 

RETURN 

END 

***** SUBROUTINE NOLIN ***** 

SUBROUTINE NONLIN (FUNC , NS IG , NE , NC , X, Y , Yl , DELTA, A, IPVT , WORK) 
IMPLICIT REAL > ' t 8(A-H,0-Z) 

DIMENSION X (NE) , Y (NE) , Yl (NE) ,A(NE,NE) , IPVT(NE) ,WORK(NE) 

EXTERNAL FUNC 

NDIM=NE 

EPSI=1.DOO/10.DOO**NSIG 

CALL NONLIO (FUNC , EPS I , NE , NC , X , DELTA, NDIM , A, Y , Y 1 , WORK , IP VT) 

RETURN 

END 


***** SUBROUTINE NOLINO ***** 

SUBROUTINE NONLIO (FUNC , EPSI , NE , NC , X, DELTA, NDIM, A, Y , Yl , WORK , IPVT) 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION X (NE) ,Y(NE) ,Y1 (NE) , IPVT (NE) ,WORK(NE) ,A(NDIM,NE) 
EXTERNAL FUNC 
NC: // OF COUNT TIMES 
DO 5 1=1, NC 
CALL FUNC (X, Y,NE) 

NE: // OF EQUATIONS 


DO 15 J= 1 , NE 

IF (DABS(Y(J)) .GT.EPSI) GO TO 25 
15 CONTINUE 
GO TO 105 
25 DO 35 J= 1 , NE 
35 Y1(J)=Y(J) 

DO 45 J= 1 , NE 

DIFF=DABS (X ( J) ) "DELTA 

IF (DABS (X (J) ) . LT . 1 . D - 12) DIFF=DELTA 

XMAM=X(J) 

X(J)=X(J) -DIFF 
CALL FUNC (X , Y , NE) 

X ( J) =XMAM 
DO 55 K= 1 , NE 

A(K, J) = ( Y 1 (K)-Y(K))/DIFF 
55 CONTINUE 
45 CONTINUE 

DO 65 J= 1 , NE 
65 Y(J)=-Y1(J) 

CALL DECOMP (NDIM,NE, A.COND, IPVT, WORK) 

CALL SOLVE (NDIM, NE, A, Y , IPVT) 

DO 75 J= 1 , NE 
X ( J) =X ( J) +Y ( J) 

75 CONTINUE 
5 CONTINUE 
105 RETURN 
END 

***** SUBROUTINE DECOMP ***** 

SUBROUTINE DECOMP (NDIM, N , A, COND , IPVT , WORK) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A (NDIM, N) ,WORK(N) , IPVT(N) 

DECOMPOSES A REAL MATRIX BY GAUSSIAN ELIMINATION, 

AND ESTIMATES THE CONDITION OF THE MATRIX. 

-COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS-, BY G. E. FORSYTHE, 
M. A. MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

USE SUBROUTINE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEM. 

INPUT. . 

NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A 
N = ORDER OF THE MATRIX 
A = MATRIX TO BE TRIANGULAR I ZED 

OUTPUT. . 

A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PREMUTED 
VERSION OF A LOWER TRIANGULAR MATRIX I~L SO THAT 
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(PERMUTATION MATRIX) *A=L*U 


COND = AN ESTIMATE OF THE CONDITION OF A. 

FOR THE LINEAR SYSTEM A*X = B , CHANGES IN A AND B 
MAY CAUSE CHANGES COND TIMES AS LARGE IN X. 

IF COND+l.O .EQ. COND , A IS SINGULAR TO WORKING 
PRECISION. COND IS SET TO 1.0D+32 IF EXACT 
SINGULARITY IS DETECTED. 

IPVT = THE PIVOT VECTOR 

IPVT(K) = THE INDEX OF THE K-TH PIVOT ROW 
IPVT (N) = (-1)** (NUMBER OF INTERCHANGES) 

WORK SPACE.. THE VECTOR WORK MUST BE DECLARED AND INCLUDED 
IN THE CALL. ITS INPUT CONTENTS ARE IGNORED. 

ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. 

THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY 
DET(A) = IPVT (N) * A (1 , 1) * A(2,2) * ... * A(N,N) 

IPVT (N) = 1 

IF (N.EQ.l) GO TO 150 
NM1=N-1 

COMPUTE THE 1-NORM OF A . 

ANORM=0 . DO 
DO 20 J=1 , N 
T=0.D0 
DO 10 1=1, N 
10 T=T+DABS (A (I , j) ) 

IF (T.GT.ANORM) ANORM=T 
20 CONTINUE 

DO GAUSSIAN ELIMINATION WITH PARTIAL 
PIVOTING. 

DO 70 K=1 , NM1 
KP1=K+1 

FIND THE PIVOT. 

M=K 

DO 30 I=KP1,N 

IF (DABS (A(l , K) ) . GT.DABS (A (M,K) ) ) M=I 
30 CONTINUE 
IPVT (K) =M 

IF (M.NE.K) IPVT (N)=- IPVT (N) 

T=A(M,K) 

A (M , K) =A(K,K) 

A(K,K)=T 

SKIP THE ELIMINATION STEP IF PIVOT IS ZERO 
IF (T.EQ.O.DO) GO TO 70 

COMPUTE THE MULTIPLIERS. 

DO 40 I=KP1 , N 
40 A ( I , K) =- A ( I , K) /T 


DO 60 J=KP1 , N 


INTERCHANGE AND ELIMINATE BY COLUMNS. 



T=A (M, J) 

A(M, J)=A(K, J) 

A(K, J)=T 

IF (T.EQ.O.DO) GO TO 60 
DO 50 I=KP 1 , N 

50 A(I, J)=A(I, J)+A(I,K)*T 
60 CONTINUE 
70 CONTINUE 


* 

* COND = (1-NORM OF A) * (AN ESTIMATE OF THE 1-NORM OF A- INVERSE) 

* THE ESTIMATE IS OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE 

* SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS 

* OF EQUATIONS, (A- TRANSPOSE) *Y = E AND A*Z = Y WHERE E 

* IS A VECTOR OF +1 OR -1 COMPONENTS CHOSEN TO CAUSS GROWTH IN Y. 

* ESTIMATE = (1-NORM OF Z) / (1-NORM OF Y) 

* 

* SOLVE (A- TRANSPOSE) *Y = E . 

DO 100 K= 1 , N 

T=0. DO 

IF (K.EQ. 1) GO TO 90 
KM1=K-1 
DO 80 1=1, KM1 
80 T=T+A(l,K)*WORK(l) 

90 EK= 1 . DO 

IF (T.LT.O.DO) EK=- 1 . DO 

IF (A(K,K) .EQ.O.DO) GO TO 160 
Al 1=A (1 , 1) 

WORK (K) =- (EK+T) /A ( 1 , 1 ) 

100 CONTINUE 

DO 120 KB=1 , NM1 
K=N-KB 
T=0.D0 
KPl=K+l 

DO 110 I=KP 1 , N 
110 t=t+a(i ,K) *work(k) 

WORK (K) =T 
M=IPVT(K) 

IF (M.EQ.K) GO TO 120 
T=WORK(M) 

WORK (M) =WORK (K) 

WORK (K) =T 
120 CONTINUE 

* 

YNORM=0 . DO 
DO 130 1=1, N 

130 YNORM=YNORM+DABS (WORK (i) ) 

Vc 

* SOLVE A*Z = Y 
CALL SOLVE (NDIM,N, A, WORK, TPVT) 

* 
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ZNORM=0 . DO 
DO 140 1 = 1, N 

140 ZNORM=ZNORM+DABS (WORK (i) ) 


* 

Vc 


i c 

150 


* 


ESTIMATE THE CONDITION. 

COND=ANORM*ZNORM/YNORM 
IF (COND.LT. l.DO) COND=1.DO 
RETURN 

1-BY-l CASE.. 

COND= 1 . DO 

IF (A (1 , 1) .NE.O.DO) RETURN 


* 

160 C0ND=1.0D32 
RETURN 
END 

SUBROUTINE SOLVE 

* 


EXACT SINGULARITY 


(NDIM,N,A,B, I PVT) 


IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION A(NDIM,N) ,B(N) ,IPVT (N) 


w SOLVES A LINEAR SYSTEM, A*X = B 

* 00 N °T SOLVE THE SYSTEM IF DECOMP HAS DETECTED SINGULARITY. 


* -COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS-, BY G. E. FORSYTHE, 

* M - A - MALCOLM, AND C. B. MOLER (PRENTICE-HALL, 1977) 

* INPUT.. 


,c NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A 

* N = ORDER OF MATRIX 

* A TR I ANGULAR IZED MATRIX OBTAINED FROM SUBROUTINE DECOMP 

* B = RIGHT HAND SIDE VECTOR 

* I PVT = PIVOT VECTOR OBTAINED FROM DECOMP 

* OUTPUT . . 

/V 

* B = SOLUTION VECTOR, X 

* 

DO THE FORWARD ELIMINATION. 

IF (N.EQ.l) GO TO 50 
NM1=N-1 
DO 20 K= 1 , NMl 
KP1=K+1 
M=IPVT(K) 

T=B (M) 

B (M) =B (K) 

B(K)=T 

DO 10 I=KP 1 , N 
10 B(I)=B(I)+A(I,K)*T 

20 CONTINUE 

NOW DO THE BACK SUBSTITUTION. 

DO 40 KB= 1 , NMl 
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KM1=N-KB 

K=KM1+1 

B(K)=B(K)/a(K,K) 

T=-B(K) 

DO 30 1=1, KM1 
30 B (I) =B (I)+A(I,K)*T 
40 CONTINUE 
50 B(1)=B(1)/a(1,1) 
RETURN 
END 
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